home *** CD-ROM | disk | FTP | other *** search
/ The World of Computer Software / The World of Computer Software.iso / srcuc.zip / ARRAY.C < prev    next >
C/C++ Source or Header  |  1991-12-20  |  64KB  |  2,263 lines
  1. /* -*-C-*-
  2.  
  3. $Header: /scheme/users/cph/src/microcode/RCS/array.c,v 9.44 1991/12/20 22:48:36 cph Exp $
  4.  
  5. Copyright (c) 1987-91 Massachusetts Institute of Technology
  6.  
  7. This material was developed by the Scheme project at the Massachusetts
  8. Institute of Technology, Department of Electrical Engineering and
  9. Computer Science.  Permission to copy this software, to redistribute
  10. it, and to use it for any purpose is granted, subject to the following
  11. restrictions and understandings.
  12.  
  13. 1. Any copy made of this software must include this copyright notice
  14. in full.
  15.  
  16. 2. Users of this software agree to make their best efforts (a) to
  17. return to the MIT Scheme project any improvements or extensions that
  18. they make, so that these may be included in future releases; and (b)
  19. to inform MIT of noteworthy uses of this software.
  20.  
  21. 3. All materials developed as a consequence of the use of this
  22. software shall duly acknowledge such use, in accordance with the usual
  23. standards of acknowledging credit in academic research.
  24.  
  25. 4. MIT has made no warrantee or representation that the operation of
  26. this software will be error-free, and MIT is under no obligation to
  27. provide any services, by way of maintenance, update, or otherwise.
  28.  
  29. 5. In conjunction with products arising from the use of this material,
  30. there shall be no use of the name of the Massachusetts Institute of
  31. Technology nor of any adaptation thereof in any advertising,
  32. promotional, or sales literature without prior written consent from
  33. MIT in each case. */
  34.  
  35. #include "scheme.h"
  36. #include "prims.h"
  37. #include "array.h"
  38. #include <math.h>
  39. #include <values.h>
  40. /* <values.h> contains some math constants */
  41.  
  42. /* ARRAY (as a scheme object)
  43.    is a usual array (in C) containing REAL numbers (float/double)
  44.    and tagged as a non-marked vector.
  45.  
  46.    Basic contents:
  47.    constructors, selectors, arithmetic operations,
  48.    conversion routines between C_Array, and Scheme_Vector
  49.  
  50.    see array.h for macros, NM_VECTOR, and extern */
  51.  
  52. /* mathematical constants */
  53. #ifdef PI
  54. #undef PI
  55. #endif
  56. #define PI           3.141592653589793238462643
  57. #define PI_OVER_2    1.570796326794896619231322
  58. #define TWOPI        6.283185307179586476925287
  59. #define SQRT_2          1.4142135623730950488
  60. #define ONE_OVER_SQRT_2  .7071067811865475244
  61. /* Abramowitz and Stegun p.3 */
  62.  
  63. REAL
  64. flonum_to_real (argument, arg_number)
  65.      fast SCHEME_OBJECT argument;
  66.      int arg_number;
  67. {
  68.   switch (OBJECT_TYPE (argument))
  69.     {
  70.     case TC_FIXNUM:
  71.       return ((REAL) (FIXNUM_TO_DOUBLE (argument)));
  72.  
  73.     case TC_BIG_FIXNUM:
  74.       if (! (BIGNUM_TO_DOUBLE_P (argument)))
  75.       error_bad_range_arg (arg_number);
  76.       return ((REAL) (bignum_to_double (argument)));
  77.  
  78.     case TC_BIG_FLONUM:
  79.       return ((REAL) (FLONUM_TO_DOUBLE (argument)));
  80.  
  81.     default:
  82.       error_wrong_type_arg (arg_number);
  83.       /* NOTREACHED */
  84.     }
  85. }
  86.  
  87. SCHEME_OBJECT
  88. allocate_array (length)
  89.      long length;
  90. {
  91. #if (REAL_IS_DEFINED_DOUBLE == 0)
  92.  
  93.   fast SCHEME_OBJECT result =
  94.     (allocate_non_marked_vector
  95.      (TC_NON_MARKED_VECTOR, ((length * REAL_SIZE) + 1), true));
  96.   FAST_MEMORY_SET (result, 1, length);
  97.   return (result);
  98.  
  99. #else /* (REAL_IS_DEFINED_DOUBLE != 0) */
  100.   
  101.   long n_words = (length * DOUBLE_SIZE);
  102.   ALIGN_FLOAT (Free);
  103.   Primitive_GC_If_Needed (n_words + 1);
  104.   {
  105.     SCHEME_OBJECT result = (MAKE_POINTER_OBJECT (TC_BIG_FLONUM, (Free)));
  106.     (*Free++) = (MAKE_OBJECT (TC_MANIFEST_NM_VECTOR, n_words));
  107.     Free += n_words;
  108.     return (result);
  109.   }
  110.  
  111. #endif /* (REAL_IS_DEFINED_DOUBLE != 0) */
  112. }
  113.  
  114. DEFINE_PRIMITIVE ("VECTOR->ARRAY", Prim_vector_to_array, 1, 1, 0)
  115. {
  116.   PRIMITIVE_HEADER (1);
  117.   CHECK_ARG (1, VECTOR_P);
  118.   {
  119.     SCHEME_OBJECT vector = (ARG_REF (1));
  120.     long length = (VECTOR_LENGTH (vector));
  121.     SCHEME_OBJECT result = (allocate_array (length));
  122.     fast SCHEME_OBJECT * scan_source = (& (VECTOR_REF (vector, 0)));
  123.     fast SCHEME_OBJECT * end_source = (scan_source + length);
  124.     fast REAL * scan_target = (ARRAY_CONTENTS (result));
  125.     while (scan_source < end_source)
  126.       (*scan_target++) = (flonum_to_real ((*scan_source++), 1));
  127.     PRIMITIVE_RETURN (result);
  128.   }
  129. }
  130.  
  131. DEFINE_PRIMITIVE ("ARRAY->VECTOR", Prim_array_to_vector, 1, 1, 0)
  132. {
  133.   PRIMITIVE_HEADER (1);
  134.   CHECK_ARG (1, ARRAY_P);
  135.   {
  136.     SCHEME_OBJECT array = (ARG_REF (1));
  137.     long length = (ARRAY_LENGTH (array));
  138.     fast REAL * scan_source = (ARRAY_CONTENTS (array));
  139.     fast REAL * end_source = (scan_source + length);
  140.     SCHEME_OBJECT result = (allocate_marked_vector (TC_VECTOR, length, true));
  141.     fast SCHEME_OBJECT * scan_result = (MEMORY_LOC (result, 1));
  142.     while (scan_source < end_source)
  143.       (*scan_result++) = (double_to_flonum ((double) (*scan_source++)));
  144.     PRIMITIVE_RETURN (result);
  145.   }
  146. }
  147.  
  148. DEFINE_PRIMITIVE ("ARRAY-ALLOCATE", Prim_array_allocate, 1,1, 0)
  149. {
  150.   fast REAL * scan;
  151.   long length;
  152.   SCHEME_OBJECT result;
  153.   PRIMITIVE_HEADER (1);
  154.  
  155.   length = (arg_nonnegative_integer (1));
  156.   result = (allocate_array (length));
  157.   for (scan = (ARRAY_CONTENTS (result)); --length >= 0; )
  158.     *scan++ = ((REAL) 0.0);
  159.   PRIMITIVE_RETURN (result);
  160. }
  161.  
  162. DEFINE_PRIMITIVE ("ARRAY-CONS-REALS", Prim_array_cons_reals, 3, 3, 0)
  163. {
  164.   PRIMITIVE_HEADER (3);
  165.   {
  166.     fast double from = (arg_real_number (1));
  167.     fast double dt = (arg_real_number (2));
  168.     long length = (arg_nonnegative_integer (3));
  169.     SCHEME_OBJECT result = (allocate_array (length));
  170.     fast REAL * scan_result = (ARRAY_CONTENTS (result));
  171.     fast int i;
  172.     for (i = 0; (i < length); i += 1)
  173.       {
  174.     (*scan_result++) = ((REAL) from);
  175.     from += dt;
  176.       }
  177.     PRIMITIVE_RETURN (result);
  178.   }
  179. }
  180.  
  181. DEFINE_PRIMITIVE ("ARRAY-LENGTH", Prim_array_length, 1, 1, 0)
  182. {
  183.   PRIMITIVE_HEADER (1);
  184.   CHECK_ARG (1, ARRAY_P);
  185.   PRIMITIVE_RETURN (LONG_TO_UNSIGNED_FIXNUM (ARRAY_LENGTH (ARG_REF (1))));
  186. }
  187.  
  188. DEFINE_PRIMITIVE ("ARRAY-REF", Prim_array_ref, 2, 2, 0)
  189. {
  190.   SCHEME_OBJECT array;
  191.   PRIMITIVE_HEADER (2);
  192.   CHECK_ARG (1, ARRAY_P);
  193.   array = (ARG_REF (1));
  194.   PRIMITIVE_RETURN
  195.     (double_to_flonum
  196.      ((double)
  197.       ((ARRAY_CONTENTS (array))
  198.        [arg_index_integer (2, (ARRAY_LENGTH (array)))])));
  199. }
  200.  
  201. DEFINE_PRIMITIVE ("ARRAY-SET!", Prim_array_set, 3, 3, 0)
  202. {
  203.   SCHEME_OBJECT array;
  204.   REAL * array_ptr;
  205.   double old_value, new_value;
  206.   PRIMITIVE_HEADER (3);
  207.   CHECK_ARG (1, ARRAY_P);
  208.   array = (ARG_REF (1));
  209.   array_ptr =
  210.     (& ((ARRAY_CONTENTS (array))
  211.     [arg_index_integer (2, (ARRAY_LENGTH (array)))]));
  212.   old_value = (*array_ptr);
  213.   new_value = (arg_real_number (3));
  214. #if (REAL_IS_DEFINED_DOUBLE == 0)
  215.   if ((new_value >= 0.0)
  216.       ? (new_value < ((double) FLT_MIN))
  217.       : (new_value > (0.0 - ((double) FLT_MIN))))
  218.     new_value = ((REAL) 0.0);
  219. #endif
  220.   (*array_ptr) = ((REAL) new_value);
  221.   PRIMITIVE_RETURN (double_to_flonum (old_value));
  222. }
  223.  
  224. /*____________________ file readers ___________
  225.   ascii and 2bint formats 
  226.   ______________________________________________*/
  227.  
  228. /* Reading data from files 
  229.    To read REAL numbers, use "lf" for double, "%f" for float 
  230.    */
  231. #if (REAL_IS_DEFINED_DOUBLE == 1)
  232. #define REALREAD  "%lf"
  233. #define REALREAD2 "%lf %lf"
  234. #else
  235. #define REALREAD  "%f"
  236. #define REALREAD2 "%f %f"
  237. #endif
  238.  
  239. static void
  240. C_Array_Read_Ascii_File (a, N, fp)          /* 16 ascii decimal digits */
  241.      REAL * a;
  242.      long N;
  243.      FILE * fp;
  244. { 
  245.   fast long i;
  246.   for (i = 0; (i < N); i += 1)
  247.     {
  248.       if ((fscanf (fp, REALREAD, (&(a[i])))) != 1)
  249.     { printf("Not enough values read ---\n last value a[%d] = % .16e \n", (i-1), a[i-1]);
  250.       error_external_return (); }
  251.     }
  252.   return;
  253. }
  254.  
  255. /* 2BINT FORMAT = integer stored in 2 consecutive bytes.
  256.    On many machines, "putw" and "getw" use 4 byte integers (C int)
  257.    so use "putc" "getc" as shown below.
  258.    */
  259.  
  260. static void
  261. C_Array_Read_2bint_File (a, N, fp)
  262.      REAL * a;
  263.      long N;
  264.      FILE * fp;
  265. {
  266.   fast long i;
  267.   fast int msd;
  268.   for (i = 0; (i < N); i += 1)
  269.     {
  270.       if (feof (fp))
  271.     error_external_return ();
  272.       msd = (getc (fp));
  273.       (a [i]) = ((REAL) ((msd << 8) | (getc (fp))));
  274.     }
  275.   return;
  276. }
  277.  
  278. DEFINE_PRIMITIVE ("ARRAY-READ-FROM-FILE", Prim_array_read_from_file, 3,3, 0)
  279. {
  280.   PRIMITIVE_HEADER (3);
  281.   CHECK_ARG (1, STRING_P);    /* 1 = filename */
  282.   /*                               2 = length of data */
  283.   CHECK_ARG (3, FIXNUM_P);    /* 3 = format of data   0=ascii 1=2bint  */
  284.   {
  285.     fast long length = (arg_nonnegative_integer (2));
  286.     fast SCHEME_OBJECT result = (allocate_array (length));
  287.     int format;
  288.     fast FILE * fp;
  289.     if ( (fp = fopen((STRING_ARG (1)), "r")) == NULL)
  290.       error_bad_range_arg (1);
  291.     
  292.     format = arg_nonnegative_integer(3);
  293.     if (format==0)
  294.       C_Array_Read_Ascii_File ((ARRAY_CONTENTS (result)), length, fp);
  295.     else if (format==1)
  296.       C_Array_Read_2bint_File ((ARRAY_CONTENTS (result)), length, fp);
  297.     else
  298.       error_bad_range_arg(3);    /* illegal format code */
  299.     
  300.     if ((fclose (fp)) != 0)
  301.       error_external_return ();
  302.     PRIMITIVE_RETURN (result);
  303.   }
  304. }
  305.  
  306. static void
  307. C_Array_Write_Ascii_File (a, N, fp)           /* 16 ascii decimal digits */
  308.      REAL * a;
  309.      long N;
  310.      FILE * fp;
  311. {
  312.   fast long i;
  313.   for (i = 0; (i < N); i += 1)
  314.     {
  315.       if (feof (fp))
  316.     error_external_return ();
  317.       fprintf (fp, "% .16e \n", a[i]);
  318.     }
  319.   return;
  320. }
  321.  
  322. DEFINE_PRIMITIVE ("ARRAY-WRITE-ASCII-FILE", Prim_array_write_ascii_file, 2, 2, 0)
  323. {
  324.   PRIMITIVE_HEADER (2);
  325.   CHECK_ARG (1, ARRAY_P);
  326.   CHECK_ARG (2, STRING_P);
  327.   {
  328.     fast SCHEME_OBJECT array = (ARG_REF (1));
  329.     fast FILE * fp = (fopen((STRING_ARG (2)), "w"));
  330.     if (fp == ((FILE *) 0))
  331.       error_bad_range_arg (2);
  332.     C_Array_Write_Ascii_File
  333.       ((ARRAY_CONTENTS (array)),
  334.        (ARRAY_LENGTH (array)),
  335.        fp);
  336.     if ((fclose (fp)) != 0)
  337.       error_external_return ();
  338.   }
  339.   PRIMITIVE_RETURN (UNSPECIFIC);
  340. }
  341.  
  342.  
  343.  
  344.  
  345. DEFINE_PRIMITIVE ("SUBARRAY-COPY!", Prim_subarray_copy, 5, 5, 0)
  346. {
  347.   PRIMITIVE_HEADER (5);
  348.   CHECK_ARG (1, ARRAY_P);    /* source array */
  349.   CHECK_ARG (2, ARRAY_P);    /* destination array */
  350.   {
  351.     REAL * source = (ARRAY_CONTENTS (ARG_REF (1)));
  352.     REAL * target = (ARRAY_CONTENTS (ARG_REF (2)));
  353.     long start_source = (arg_nonnegative_integer (3));
  354.     long start_target = (arg_nonnegative_integer (4));
  355.     long n_elements = (arg_nonnegative_integer (5));
  356.     if ((start_source + n_elements) > (ARRAY_LENGTH (ARG_REF (1))))
  357.       error_bad_range_arg (3);
  358.     if ((start_target + n_elements) > (ARRAY_LENGTH (ARG_REF (2))))
  359.       error_bad_range_arg (4);
  360.     {
  361.       fast REAL * scan_source = (source + start_source);
  362.       fast REAL * end_source = (scan_source + n_elements);
  363.       fast REAL * scan_target = (target + start_target);
  364.       while (scan_source < end_source)
  365.     (*scan_target++) = (*scan_source++);
  366.     }
  367.   }
  368.   PRIMITIVE_RETURN (UNSPECIFIC);
  369. }
  370.  
  371. DEFINE_PRIMITIVE ("ARRAY-REVERSE!", Prim_array_reverse, 1, 1, 0)
  372. {
  373.   PRIMITIVE_HEADER (1);
  374.   CHECK_ARG (1, ARRAY_P);
  375.   {
  376.     SCHEME_OBJECT array = (ARG_REF (1));
  377.     long length = (ARRAY_LENGTH (array));
  378.     long half_length = (length / 2);
  379.     fast REAL * array_ptr = (ARRAY_CONTENTS (array));
  380.     fast long i;
  381.     fast long j;
  382.     for (i = 0, j = (length - 1); (i < half_length); i += 1, j -= 1)
  383.       {
  384.     fast REAL Temp = (array_ptr [j]);
  385.     (array_ptr [j]) = (array_ptr [i]);
  386.     (array_ptr [i]) = Temp;
  387.       }
  388.   }
  389.   PRIMITIVE_RETURN (UNSPECIFIC);
  390. }
  391.  
  392. DEFINE_PRIMITIVE ("ARRAY-TIME-REVERSE!", Prim_array_time_reverse, 1, 1, 0)
  393. {
  394.   void C_Array_Time_Reverse ();
  395.   PRIMITIVE_HEADER (1);
  396.   CHECK_ARG (1, ARRAY_P);
  397.   C_Array_Time_Reverse
  398.     ((ARRAY_CONTENTS (ARG_REF (1))), (ARRAY_LENGTH (ARG_REF (1))));
  399.   PRIMITIVE_RETURN (UNSPECIFIC);
  400. }
  401.  
  402. /* time-reverse
  403.    x[0] remains fixed. (time=0)
  404.    x[i] swapped with x[n-i]. (mirror image around x[0]) */
  405.  
  406. void
  407. C_Array_Time_Reverse (x,n)
  408.      REAL *x;
  409.      long n;
  410. { long i, ni, n2;
  411.   REAL xt;
  412.   if ((n % 2) == 0)        /* even length */
  413.   { n2 = (n/2);
  414.     for (i=1; i<n2; i++)    /* i=1,2,..,n/2-1 */
  415.     {  ni = n-i;
  416.        xt    = x[i];
  417.        x[i]  = x[ni];
  418.        x[ni] = xt; }}
  419.   else                /* odd length */
  420.   { n2 = (n+1)/2;        /* (n+1)/2 = (n-1)/2 + 1 */
  421.     for (i=1; i<n2; i++)    /* i=1,2,..,(n-1)/2 */
  422.     {  ni = n-i;
  423.        xt   = x[i];
  424.        x[i] = x[ni];
  425.        x[ni] = xt; }}
  426. }
  427.  
  428. /* The following is smart
  429.    and avoids computation when offset or scale are degenerate 0,1 */
  430.  
  431. DEFINE_PRIMITIVE ("SUBARRAY-OFFSET-SCALE!", Prim_subarray_offset_scale, 5, 5, 0)
  432. {
  433.   long i, at, m,mplus;
  434.   REAL *a, offset,scale;
  435.   PRIMITIVE_HEADER (5);
  436.   CHECK_ARG (1, ARRAY_P);
  437.   CHECK_ARG (2, FIXNUM_P);
  438.   CHECK_ARG (3, FIXNUM_P);
  439.   a = ARRAY_CONTENTS(ARG_REF(1));
  440.   at = arg_nonnegative_integer(2); /*       at = starting index             */
  441.   m  = arg_nonnegative_integer(3); /*       m  = number of points to change */
  442.   mplus = at + m;
  443.   if (mplus > (ARRAY_LENGTH(ARG_REF(1)))) error_bad_range_arg(3);
  444.   offset = (arg_real (4));
  445.   scale = (arg_real (5));
  446.   if ((offset == 0.0) && (scale == 1.0))
  447.     ;                /* be smart */
  448.   else if (scale == 0.0)
  449.     for (i=at; i<mplus; i++)  a[i] = offset;
  450.   else if (offset == 0.0)
  451.     for (i=at; i<mplus; i++)  a[i] = scale * a[i];
  452.   else if (scale == 1.0)
  453.     for (i=at; i<mplus; i++)  a[i] = offset + a[i];
  454.   else
  455.     for (i=at; i<mplus; i++)  a[i] = offset + scale * a[i];
  456.   PRIMITIVE_RETURN (UNSPECIFIC);
  457. }
  458.  
  459. DEFINE_PRIMITIVE ("COMPLEX-SUBARRAY-COMPLEX-SCALE!", Prim_complex_subarray_complex_scale, 6,6, 0)
  460. { long i, at,m,mplus;
  461.   REAL *a,*b;            /* (a,b) = (real,imag) arrays */
  462.   double temp, minus_y,  x, y;    /* (x,y) = (real,imag) scale */
  463.   PRIMITIVE_HEADER (6);
  464.   CHECK_ARG (1, ARRAY_P);
  465.   CHECK_ARG (2, ARRAY_P);
  466.   at = (arg_nonnegative_integer (3)); /* starting index */
  467.   m  = (arg_nonnegative_integer (4)); /* number of points to change */
  468.   mplus = at + m;
  469.   if (mplus > (ARRAY_LENGTH(ARG_REF(1)))) error_bad_range_arg(4);
  470.   x = (arg_real_number (5));
  471.   y = (arg_real_number (6));
  472.   a = ARRAY_CONTENTS(ARG_REF(1));
  473.   b = ARRAY_CONTENTS(ARG_REF(2));
  474.   if ((ARRAY_LENGTH(ARG_REF(1))) != (ARRAY_LENGTH(ARG_REF(2))))
  475.     error_bad_range_arg(2);
  476.   if (x==0.0)            /* imaginary only */
  477.     if       (y==0.0)
  478.       for (i=at; i<mplus; i++)
  479.       { a[i] = 0.0;
  480.     b[i] = 0.0; }
  481.     else if  (y==1.0)
  482.       for (i=at; i<mplus; i++)
  483.       { temp = b[i];
  484.     b[i] = a[i];
  485.     a[i] = (-temp); }
  486.     else if  (y==-1.0)
  487.       for (i=at; i<mplus; i++)
  488.       { temp = b[i];
  489.     b[i] = (-a[i]);
  490.     a[i] = temp; }
  491.     else
  492.     { minus_y = (-y);
  493.       for (i=at; i<mplus; i++)
  494.       { temp =               y * ((double) a[i]);
  495.     a[i] = (REAL) (minus_y * ((double) b[i]));
  496.     b[i] = (REAL) temp; }}
  497.   else if (y==0.0)        /* real only */
  498.     if (x==1.0) ;
  499.     else for (i=at; i<mplus; i++)
  500.     { a[i] = (REAL) (x * ((double) a[i]));
  501.       b[i] = (REAL) (x * ((double) b[i])); }
  502.   else                /* full complex scale */
  503.     for (i=at; i<mplus; i++)
  504.     { temp =         ((double) a[i])*x - ((double) b[i])*y;
  505.       b[i] = (REAL) (((double) b[i])*x + ((double) a[i])*y);
  506.       a[i] = (REAL) temp; }
  507.   PRIMITIVE_RETURN (UNSPECIFIC);
  508. }
  509.  
  510. /* Accumulate
  511.    using combinators              *
  512.    corresponding type codes       1 */
  513.  
  514. DEFINE_PRIMITIVE ("COMPLEX-SUBARRAY-ACCUMULATE!", Prim_complex_subarray_accumulate, 6,6, 0)
  515. { long  at,m,mplus, tc, i;
  516.   REAL *a,*b;            /* (a,b) = (real,imag) input arrays */
  517.   REAL *c;    /* result = output array of length 2, holds a complex number */
  518.   double x, y, temp;
  519.   PRIMITIVE_HEADER (6);
  520.   CHECK_ARG (1, ARRAY_P);    /* a = input array (real) */
  521.   CHECK_ARG (2, ARRAY_P);    /* b = input array (imag) */
  522.   a = ARRAY_CONTENTS(ARG_REF(1));
  523.   b = ARRAY_CONTENTS(ARG_REF(2));
  524.   if ((ARRAY_LENGTH(ARG_REF(1))) != (ARRAY_LENGTH(ARG_REF(2))))
  525.     error_bad_range_arg(2);
  526.   tc = arg_nonnegative_integer(3); /*       tc = type code 0 or 1            */
  527.   at = arg_nonnegative_integer(4); /*       at = starting index              */
  528.   m  = arg_nonnegative_integer(5); /*       m  = number of points to process */
  529.   CHECK_ARG (6, ARRAY_P);    /* c = output array of length 2 */
  530.   c = ARRAY_CONTENTS(ARG_REF(6));
  531.   if ((ARRAY_LENGTH(ARG_REF(6))) != 2) error_bad_range_arg(6);
  532.   mplus = at + m;
  533.   if (mplus > (ARRAY_LENGTH(ARG_REF(1)))) error_bad_range_arg(5);
  534.   if (tc==1)
  535.   { x = 1.0;            /* real part of accumulator */
  536.     y = 0.0;            /* imag part of accumulator */
  537.     for (i=at;i<mplus;i++)
  538.     { temp = ((double) a[i])*x - ((double) b[i])*y;
  539.       y    = ((double) b[i])*x + ((double) a[i])*y;
  540.       x    = temp; }
  541.   }
  542.   else
  543.     error_bad_range_arg(3);
  544.   c[0] = ((REAL) x);        /* mechanism for returning complex number */
  545.   c[1] = ((REAL) y);        /* do not use lists, avoid heap pointer   */
  546.   PRIMITIVE_RETURN (UNSPECIFIC);
  547. }
  548.  
  549. DEFINE_PRIMITIVE ("CS-ARRAY-TO-COMPLEX-ARRAY!", Prim_cs_array_to_complex_array, 3, 3, 0)
  550. { long n,n2,n2_1, i;
  551.   REAL *a, *b,*c;
  552.   PRIMITIVE_HEADER (3);
  553.   CHECK_ARG (1, ARRAY_P);
  554.   CHECK_ARG (2, ARRAY_P);
  555.   CHECK_ARG (3, ARRAY_P);
  556.   a = ARRAY_CONTENTS(ARG_REF(1));
  557.   n = ARRAY_LENGTH(ARG_REF(1));
  558.   b = ARRAY_CONTENTS(ARG_REF(2));
  559.   c = ARRAY_CONTENTS(ARG_REF(3));
  560.   if (n!=(ARRAY_LENGTH(ARG_REF(2)))) error_bad_range_arg(2);
  561.   if (n!=(ARRAY_LENGTH(ARG_REF(3)))) error_bad_range_arg(3);
  562.   b[0] = a[0];
  563.   c[0] = 0.0;
  564.   n2 = n/2;            /* integer division truncates down */
  565.   n2_1 = n2+1;
  566.   if (2*n2 == n)        /* even length, n2 is only real */
  567.   { b[n2]  = a[n2];  c[n2]  = 0.0; }
  568.   else /* odd length, make the loop include the n2 index */
  569.   { n2   = n2+1;
  570.     n2_1 = n2; }
  571.   for (i=1; i<n2; i++)   { b[i] = a[i];
  572.                c[i] = a[n-i]; }
  573.   for (i=n2_1; i<n; i++) { b[i] =  a[n-i];
  574.                c[i] = (-a[i]); }
  575.   PRIMITIVE_RETURN (UNSPECIFIC);
  576. }
  577.  
  578. DEFINE_PRIMITIVE ("CS-ARRAY-MULTIPLY-INTO-SECOND-ONE!", Prim_cs_array_multiply_into_second_one, 2, 2, 0)
  579. { long n,n2;
  580.   REAL *a, *b;
  581.   void cs_array_multiply_into_second_one();
  582.   PRIMITIVE_HEADER (2);
  583.   CHECK_ARG (1, ARRAY_P);
  584.   CHECK_ARG (2, ARRAY_P);
  585.   a = ARRAY_CONTENTS(ARG_REF(1));
  586.   n = ARRAY_LENGTH(ARG_REF(1));
  587.   b = ARRAY_CONTENTS(ARG_REF(2));
  588.   if (n!=(ARRAY_LENGTH(ARG_REF(2)))) error_bad_range_arg(2);
  589.   n2 = n/2;            /* integer division truncates down */
  590.   cs_array_multiply_into_second_one(a,b, n,n2);
  591.   PRIMITIVE_RETURN (UNSPECIFIC);
  592. }
  593.  
  594. void
  595. cs_array_multiply_into_second_one (a,b, n,n2)
  596.      REAL *a, *b;
  597.      long n,n2;
  598. {
  599.   REAL temp;
  600.   long i,ni;
  601.   b[0]   = a[0]  * b[0];
  602.   if (2*n2 == n)        /* even length, n2 is only real */
  603.     b[n2]  = a[n2] * b[n2];
  604.   else
  605.     n2 = n2+1; /* odd length, make the loop include the n2 index */
  606.   for (i=1; i<n2; i++)
  607.     {
  608.       ni = n-i;
  609.       temp   = a[i]*b[i]   -  a[ni]*b[ni]; /* real part */
  610.       b[ni]  = a[i]*b[ni]  +  a[ni]*b[i]; /*  imag part */
  611.       b[i]   = temp;
  612.     }
  613. }
  614.  
  615. DEFINE_PRIMITIVE ("CS-ARRAY-DIVIDE-INTO-XXX!", Prim_cs_array_divide_into_xxx, 4, 4, 0)
  616. {
  617.   long n,n2, one_or_two;
  618.   REAL *a, *b, inf;
  619.   void cs_array_divide_into_z();
  620.   PRIMITIVE_HEADER (4);
  621.   CHECK_ARG (1, ARRAY_P);
  622.   CHECK_ARG (2, ARRAY_P);
  623.   inf = (arg_real (3));
  624.   /* where to store result of division */
  625.   one_or_two = (arg_nonnegative_integer (4));
  626.   a = ARRAY_CONTENTS(ARG_REF(1));
  627.   b = ARRAY_CONTENTS(ARG_REF(2));
  628.   n = ARRAY_LENGTH(ARG_REF(1));
  629.   if (n!=(ARRAY_LENGTH(ARG_REF(2)))) error_bad_range_arg(2);
  630.   n2 = n/2;            /* integer division truncates down */
  631.   if (one_or_two == 1)
  632.     cs_array_divide_into_z(a,b, a,  n,n2, inf);
  633.   else if (one_or_two == 2)
  634.     cs_array_divide_into_z(a,b, b,  n,n2, inf);
  635.   else
  636.     error_bad_range_arg(4);
  637.   PRIMITIVE_RETURN (UNSPECIFIC);
  638. }
  639.  
  640. void
  641. cs_array_divide_into_second_one (a,b, n,n2,inf)    /* used in image.c */
  642.      REAL *a,*b, inf;
  643.      long n,n2;
  644. {
  645.   void cs_array_divide_into_z ();
  646.   cs_array_divide_into_z (a,b, b, n,n2,inf);
  647. }
  648.  
  649. void
  650. cs_array_divide_into_z (a,b, z, n,n2, inf) /* z can be either a or b */
  651.      REAL *a,*b,*z, inf;
  652.      long n,n2;
  653. { long i,ni;
  654.   REAL temp, radius;
  655.  
  656.   if (b[0] == 0.0)
  657.     if (a[0] == 0.0) z[0] = 1.0;
  658.     else             z[0] = a[0] * inf;
  659.   else               z[0] = a[0] / b[0];
  660.  
  661.   if (2*n2 == n)        /* even length, n2 is only real */
  662.     if (b[n2] == 0.0)
  663.       if (a[n2] == 0.0) z[n2] = 1.0;
  664.       else              z[n2] = a[n2] * inf;
  665.     else                z[n2] = a[n2] / b[n2];
  666.   else
  667.     n2 = n2+1; /* odd length, make the loop include the n2 index */
  668.  
  669.   for (i=1; i<n2; i++)
  670.   { ni = n-i;
  671.     radius = b[i]*b[i] + b[ni]*b[ni]; /* b^2 denominator = real^2 + imag^2 */
  672.  
  673.     if (radius == 0.0) {
  674.       if (a[i]  == 0.0) z[i]  = 1.0;
  675.       else              z[i]  = a[i] * inf;
  676.       if (a[ni] == 0.0) z[ni] = 1.0;
  677.       else              z[ni] = a[ni] * inf; }
  678.     else {
  679.       temp  = a[i]*b[i]    +  a[ni]*b[ni];
  680.       z[ni] = (a[ni]*b[i]  -  a[i]*b[ni]) / radius; /* imag part */
  681.       z[i]  = temp                        / radius; /* real part */
  682.     }}
  683. }
  684.  
  685. /* ARRAY-UNARY-FUNCTION!
  686.    apply unary-function elementwise on array
  687.    Available functions : */
  688.  
  689. void
  690. REALabs (a,b)
  691.      REAL *a,*b;
  692. {
  693.   (*b) = ( (REAL) fabs( (double) (*a)) );
  694. }
  695.  
  696. void
  697. REALexp (a,b)
  698.      REAL *a,*b;
  699. {
  700.   fast double y;
  701.   if ((y = exp((double) (*a))) == HUGE)
  702.     error_bad_range_arg (1);    /* OVERFLOW */
  703.   (*b) = ((REAL) y);
  704. }
  705.  
  706. void
  707. REALlog (a,b)
  708.      REAL *a,*b;
  709. {
  710.   if ((*a) < 0.0)
  711.     error_bad_range_arg(1);    /* log(negative) */
  712.   (*b) = ( (REAL) log( (double) (*a)) );
  713. }
  714.  
  715. void
  716. REALtruncate (a,b)
  717.      REAL *a,*b;        /* towards zero */
  718. {
  719.   double integral_part, modf();
  720.   modf( ((double) (*a)), &integral_part);
  721.   (*b) = ( (REAL) integral_part);
  722. }
  723.  
  724. void
  725. REALround (a,b)
  726.      REAL *a,*b;        /* towards nearest integer */
  727. {
  728.   double integral_part, modf();
  729.   if ((*a) >= 0.0)        /* It may be faster to look at the sign
  730.                    of mantissa, and dispatch */
  731.     modf( ((double) ((*a)+0.5)), &integral_part);
  732.   else
  733.     modf( ((double) ((*a)-0.5)), &integral_part);
  734.   (*b) = ( (REAL) integral_part);
  735. }
  736.  
  737. void
  738. REALsquare (a,b)
  739.      REAL *a,*b;
  740. {
  741.   (*b) = ( (REAL) ((*a) * (*a)) );
  742. }
  743.  
  744. void
  745. REALsqrt (a,b)
  746.      REAL *a,*b;
  747. {
  748.   if ((*a) < 0.0)
  749.     error_bad_range_arg(1);    /* sqrt(negative) */
  750.   (*b) = ( (REAL) sqrt( (double) (*a)) );
  751. }
  752.  
  753. void
  754. REALsin (a,b)
  755.      REAL *a,*b;
  756. {
  757.   (*b) = ( (REAL) sin( (double) (*a)) );
  758. }
  759.  
  760. void
  761. REALcos (a,b)
  762.      REAL *a,*b;
  763. {
  764.   (*b) = ( (REAL) cos( (double) (*a)) );
  765. }
  766.  
  767. void
  768. REALtan (a,b)
  769.      REAL *a,*b;
  770. {
  771.   (*b) = ( (REAL) tan( (double) (*a)) );
  772. }
  773.  
  774. void
  775. REALasin (a,b)
  776.      REAL *a,*b;
  777. {
  778.   (*b) = ( (REAL) asin( (double) (*a)) );
  779. }
  780.  
  781. void
  782. REALacos (a,b)
  783.      REAL *a,*b;
  784. {
  785.   (*b) = ( (REAL) acos( (double) (*a)) );
  786. }
  787.  
  788. void
  789. REALatan (a,b)
  790.      REAL *a,*b;
  791. {
  792.   (*b) = ( (REAL) atan( (double) (*a)) );
  793. }
  794.  
  795. void
  796. REALgamma (a,b)
  797.      REAL *a,*b;
  798. {
  799.   fast double y;
  800.   if ((y = gamma(((double) (*a)))) > LN_MAXDOUBLE)
  801.     error_bad_range_arg(1);    /* gamma( non-positive integer ) */
  802.   (*b) = ((REAL) (signgam * exp(y))); /* see HPUX Section 3 */
  803. }
  804.  
  805. void
  806. REALerf (a,b)
  807.      REAL *a,*b;
  808. {
  809.   (*b) = ( (REAL) erf((double) (*a)) );
  810. }
  811.  
  812. void
  813. REALerfc (a,b)
  814.      REAL *a,*b;
  815. {
  816.   (*b) = ( (REAL) erfc((double) (*a)) );
  817. }
  818.  
  819. void
  820. REALbessel1 (order,a,b)        /* Bessel of first kind */
  821.      long order;
  822.      REAL *a,*b;
  823. {
  824.   if (order == 0)
  825.     (*b) = ( (REAL) j0((double) (*a)) );
  826.   if (order == 1)
  827.     (*b) = ( (REAL) j1((double) (*a)) );
  828.   else
  829.     (*b) = ( (REAL) jn(((int) order), ((double) (*a))) );
  830. }
  831.  
  832. void
  833. REALbessel2 (order,a,b)        /* Bessel of second kind */
  834.      long order;
  835.      REAL *a,*b;
  836. {
  837.   if ((*a) <= 0.0)
  838.     error_bad_range_arg(1);    /* Blows Up */
  839.   if (order == 0)
  840.     (*b) = ( (REAL) y0((double) (*a)) );
  841.   if (order == 1)
  842.     (*b) = ( (REAL) y1((double) (*a)) );
  843.   else
  844.     (*b) = ( (REAL) yn(((int) order), ((double) (*a))) );
  845. }
  846.  
  847. /* Table to store the available unary-functions.
  848.    Also some binary functions at the end -- not available right now.
  849.    The (1 and 2)s denote the numofargs (1 for unary 2 for binary) */
  850.  
  851. struct array_func_table
  852. {
  853.   long numofargs;
  854.   void (*func)();
  855. }
  856. Array_Function_Table [] =
  857. {
  858.   1, REALabs,            /*0*/
  859.   1, REALexp,            /*1*/
  860.   1, REALlog,            /*2*/
  861.   1, REALtruncate,        /*3*/
  862.   1, REALround,            /*4*/
  863.   1, REALsquare,        /*5*/
  864.   1, REALsqrt,            /*6*/
  865.   1, REALsin,            /*7*/
  866.   1, REALcos,            /*8*/
  867.   1, REALtan,            /*9*/
  868.   1, REALasin,            /*10*/
  869.   1, REALacos,            /*11*/
  870.   1, REALatan,            /*12*/
  871.   1, REALgamma,            /*13*/
  872.   1, REALerf,            /*14*/
  873.   1, REALerfc,            /*15*/
  874.   2, REALbessel1,        /*16*/
  875.   2, REALbessel2        /*17*/
  876.   };
  877.  
  878. #define MAX_ARRAY_FUNCTC 17
  879.  
  880. /* array-unary-function!    could be called        array-operation-1!
  881.    following the naming convention for other similar procedures
  882.    but it is specialized to mappings only, so we have special name. */
  883.  
  884. DEFINE_PRIMITIVE ("ARRAY-UNARY-FUNCTION!", Prim_array_unary_function, 2,2, 0)
  885. {
  886.   long n, i;
  887.   REAL *a,*b;
  888.   long tc;
  889.   void (*f)();
  890.   PRIMITIVE_HEADER (2);
  891.   CHECK_ARG (1, ARRAY_P);    /*      a  = input (and output) array    */
  892.   CHECK_ARG (2, FIXNUM_P);    /*      tc = type code                   */
  893.   tc = arg_nonnegative_integer(2);
  894.   if (tc > MAX_ARRAY_FUNCTC) error_bad_range_arg(2);
  895.   f = ((Array_Function_Table[tc]).func);
  896.   if (1 != (Array_Function_Table[tc]).numofargs) error_wrong_type_arg(2);
  897.   /* check it is a unary function */
  898.   a = ARRAY_CONTENTS(ARG_REF(1));
  899.   b = a;
  900.   n = ARRAY_LENGTH(ARG_REF(1));
  901.   for (i=0; i<n; i++)
  902.     (*f) ( &(a[i]), &(b[i]) );    /* a into b */
  903.   PRIMITIVE_RETURN (UNSPECIFIC);
  904. }
  905.  
  906. /* Accumulate
  907.    using combinators              +  or  *
  908.    corresponding type codes       0      1 */
  909.  
  910. DEFINE_PRIMITIVE ("SUBARRAY-ACCUMULATE", Prim_subarray_accumulate, 4,4, 0)
  911. {
  912.   long at,m,mplus, tc, i;
  913.   REAL *a;
  914.   double result;
  915.   PRIMITIVE_HEADER (4);
  916.   CHECK_ARG (1, ARRAY_P);    /*           a = input array                 */
  917.   a  = ARRAY_CONTENTS(ARG_REF(1));
  918.   tc = arg_nonnegative_integer(2); /*       tc = type code 0 or 1            */
  919.   at = arg_nonnegative_integer(3); /*       at = starting index              */
  920.   m  = arg_nonnegative_integer(4); /*       m  = number of points to process */
  921.   mplus = at + m;
  922.   if (mplus > (ARRAY_LENGTH(ARG_REF(1)))) error_bad_range_arg(4);
  923.   if (tc == 0)
  924.     {
  925.       result = 0.0;
  926.       for (i=at;i<mplus;i++)
  927.     result = result + ((double) a[i]);
  928.     }
  929.   else if (tc == 1)
  930.     {
  931.       result = 1.0;
  932.       for (i=at;i<mplus;i++)
  933.     result = result * ((double) a[i]);
  934.     }
  935.   else
  936.     error_bad_range_arg (2);
  937.   PRIMITIVE_RETURN (double_to_flonum (result));
  938. }
  939.  
  940. /* The following searches for value within tolerance
  941.    starting from index=from in array.
  942.    Returns first index where match occurs (useful for finding zeros). */
  943.  
  944. DEFINE_PRIMITIVE ("ARRAY-SEARCH-VALUE-TOLERANCE-FROM", Prim_array_search_value_tolerance_from, 4, 4, 0)
  945. {
  946.   SCHEME_OBJECT array;
  947.   fast long length;
  948.   fast REAL * a;
  949.   fast REAL value;        /* value to search for */
  950.   fast double tolerance;    /* tolerance allowed */
  951.   PRIMITIVE_HEADER (4);
  952.   CHECK_ARG (1, ARRAY_P);
  953.   array = (ARG_REF (1));
  954.   length = (ARRAY_LENGTH (array));
  955.   a = (ARRAY_CONTENTS (array));
  956.   value = (arg_real (2));
  957.   tolerance = (arg_real (3));
  958.   {
  959.     fast long i;
  960.     for (i = (arg_index_integer (4, length)); i<length; i+=1)
  961.       if (tolerance >= (fabs ((double) ((a [i]) - value))))
  962.     PRIMITIVE_RETURN (LONG_TO_UNSIGNED_FIXNUM (i));
  963.   }
  964.   PRIMITIVE_RETURN (SHARP_F);
  965. }
  966.  
  967. DEFINE_PRIMITIVE ("SUBARRAY-MIN-MAX-INDEX", Prim_subarray_min_max_index, 3, 3, 0)
  968. {
  969.   PRIMITIVE_HEADER (3);
  970.   CHECK_ARG (1, ARRAY_P);
  971.   {
  972.     REAL * a = (ARRAY_CONTENTS (ARG_REF (1)));
  973.     long at = (arg_nonnegative_integer (2)); /* starting index */
  974.     long m  = (arg_nonnegative_integer (3)); /* number of points to process */
  975.     long mplus = (at + m);
  976.     long nmin;
  977.     long nmax;
  978.     if (mplus > (ARRAY_LENGTH (ARG_REF (1))))
  979.       error_bad_range_arg (3);
  980.     
  981.     C_Array_Find_Min_Max ((& (a [at])), m, (&nmin), (&nmax));
  982.     nmin = nmin + at;        /* offset appropriately */
  983.     nmax = nmax + at;
  984.     
  985.     PRIMITIVE_RETURN
  986.       (cons ((LONG_TO_FIXNUM (nmin)),
  987.          (cons ((LONG_TO_FIXNUM (nmax)),
  988.             EMPTY_LIST))));
  989.   }
  990. }
  991.  
  992. void
  993. C_Array_Find_Min_Max (x, n, nmin, nmax)
  994.      fast REAL * x;
  995.      fast long n;
  996.      long * nmin;
  997.      long * nmax;
  998. { REAL *xold = x;
  999.   register REAL xmin, xmax;
  1000.   register long nnmin, nnmax;
  1001.   register long count;
  1002.  
  1003.   nnmin = nnmax = 0;
  1004.   xmin = xmax = *x++;
  1005.   n--;
  1006.   count = 1;
  1007.   if(n>0)
  1008.   {
  1009.     do {
  1010.       if(*x < xmin) {
  1011.     nnmin = count++ ;
  1012.     xmin = *x++ ;
  1013.       } else if(*x > xmax) {
  1014.     nnmax = count++ ;
  1015.     xmax = *x++ ;
  1016.       } else {
  1017.     count++ ;
  1018.     x++ ;
  1019.       }
  1020.     } while( --n > 0 ) ;
  1021.   }
  1022.   *nmin = nnmin ;
  1023.   *nmax = nnmax ;
  1024. }
  1025.  
  1026.  
  1027. /* array-average
  1028.    can be done with (array-reduce +) and division by array-length.
  1029.    But there is also this C primitive.
  1030.    Keep it around, may be useful someday. */
  1031.  
  1032. /* Computes the average in pieces, so as to reduce
  1033.    roundoff smearing in cumulative sum.
  1034.    example= first huge positive numbers, then small nums,
  1035.    then huge negative numbers. */
  1036.  
  1037. static void
  1038. C_Array_Find_Average (Array, Length, pAverage)
  1039.      long Length;
  1040.      REAL * Array;
  1041.      REAL * pAverage;
  1042. {
  1043.   long i;
  1044.   long array_index;
  1045.   REAL average_n, sum;
  1046.  
  1047.   average_n = 0.0;
  1048.   array_index = 0;
  1049.   while (array_index < Length)
  1050.     {
  1051.       sum = 0.0;
  1052.       for (i=0;((array_index<Length) && (i<2000));i++) {
  1053.     sum += Array[array_index];
  1054.     array_index++;
  1055.       }
  1056.       average_n += (sum / ((REAL) Length));
  1057.     }
  1058.   (*pAverage) = average_n;
  1059.   return;
  1060. }
  1061.  
  1062. DEFINE_PRIMITIVE ("ARRAY-AVERAGE", Prim_array_find_average, 1, 1, 0)
  1063. {
  1064.   SCHEME_OBJECT array;
  1065.   REAL average;
  1066.   PRIMITIVE_HEADER (1);
  1067.   CHECK_ARG (1, ARRAY_P);
  1068.   array = (ARG_REF (1));
  1069.   C_Array_Find_Average
  1070.     ((ARRAY_CONTENTS (array)),
  1071.      (ARRAY_LENGTH (array)),
  1072.      (&average));
  1073.   PRIMITIVE_RETURN (double_to_flonum ((double) average));
  1074. }
  1075.  
  1076. void
  1077. C_Array_Make_Histogram (Array, Length, Histogram, npoints)
  1078.      REAL Array[], Histogram[];
  1079.      long Length, npoints;
  1080. {
  1081.   REAL Max, Min, Offset, Scale;
  1082.   long i, nmin, nmax, index;
  1083.   C_Array_Find_Min_Max (Array, Length, (&nmin), (&nmax));
  1084.   Min = (Array [nmin]);
  1085.   Max = (Array [nmax]);
  1086.   Find_Offset_Scale_For_Linear_Map
  1087.     (Min, Max, 0.0, ((REAL) npoints), (&Offset), (&Scale));
  1088.   for (i = 0; (i < npoints); i += 1)
  1089.     (Histogram [i]) = 0.0;
  1090.   for (i = 0; (i < Length); i += 1)
  1091.     {
  1092.       /* Everything from 0 to 1 maps to bin 0, and so on */
  1093.       index = ((long) (floor ((double) ((Scale * (Array [i])) + Offset))));
  1094.       /* max that won't floor to legal array index */
  1095.       if (index == npoints)
  1096.     index = (index - 1);
  1097.       (Histogram [index]) += 1.0;
  1098.     }
  1099.   return;
  1100. }
  1101.  
  1102. DEFINE_PRIMITIVE ("ARRAY-MAKE-HISTOGRAM", Prim_array_make_histogram, 2, 2, 0)
  1103. {
  1104.   PRIMITIVE_HEADER (2);
  1105.   CHECK_ARG (1, ARRAY_P);
  1106.   {
  1107.     fast SCHEME_OBJECT array = (ARG_REF (1));
  1108.     long length = (ARRAY_LENGTH (array));
  1109.     long npoints = (arg_integer_in_range (2, 1, ((2 * length) + 1)));
  1110.     SCHEME_OBJECT result = (allocate_array (npoints));
  1111.     C_Array_Make_Histogram
  1112.       ((ARRAY_CONTENTS (array)),
  1113.        length,
  1114.        (ARRAY_CONTENTS (result)),
  1115.        npoints);
  1116.     PRIMITIVE_RETURN (result);
  1117.   }
  1118. }
  1119.  
  1120. /* The following geometrical map is slightly tricky. */
  1121. void
  1122. Find_Offset_Scale_For_Linear_Map (Min, Max, New_Min, New_Max, Offset, Scale)
  1123.      REAL Min, Max, New_Min, New_Max, *Offset, *Scale;
  1124. {
  1125.   if (Min != Max)
  1126.     {
  1127.       (*Scale)  = ((New_Max - New_Min) / (Max - Min));
  1128.       (*Offset) = (New_Min - ((*Scale) * Min));
  1129.     }
  1130.   else
  1131.     {
  1132.       (*Scale) =
  1133.     ((Max == 0.0)
  1134.      ? 0.0
  1135.      : (0.25 * (mabs ((New_Max - New_Min) / Max))));
  1136.       (*Offset) = ((New_Max + New_Min) / 2.0);
  1137.     }
  1138.   return;
  1139. }
  1140.  
  1141.  
  1142. DEFINE_PRIMITIVE ("ARRAY-CLIP-MIN-MAX!", Prim_array_clip_min_max, 3, 3, 0)
  1143. {
  1144.   PRIMITIVE_HEADER (3);
  1145.   CHECK_ARG (1, ARRAY_P);
  1146.   {
  1147.     SCHEME_OBJECT array = (ARG_REF (1));
  1148.     REAL xmin = (arg_real (2));
  1149.     REAL xmax = (arg_real (3));
  1150.     long Length = (ARRAY_LENGTH (array));
  1151.     REAL * From_Here = (ARRAY_CONTENTS (array));
  1152.     REAL * To_Here = From_Here;
  1153.     long i;
  1154.     if (xmin > xmax)
  1155.       error_bad_range_arg (3);
  1156.     for (i = 0; (i < Length); i += 1)
  1157.       {
  1158.     if ((*From_Here) < xmin)
  1159.       (*To_Here++) = xmin;
  1160.     else if ((*From_Here) > xmax)
  1161.       (*To_Here++) = xmax;
  1162.     else
  1163.       (*To_Here++) = (*From_Here);
  1164.     From_Here += 1;
  1165.       }
  1166.   }
  1167.   PRIMITIVE_RETURN (UNSPECIFIC);
  1168. }
  1169.  
  1170. /* complex-array-operation-1!
  1171.    groups together procedures that use 1 complex-array
  1172.    and store the result in place */
  1173.  
  1174. DEFINE_PRIMITIVE ("COMPLEX-ARRAY-OPERATION-1!", Prim_complex_array_operation_1, 3,3, 0)
  1175. { long n, i, opcode;
  1176.   REAL *a, *b;
  1177.   void complex_array_to_polar(), complex_array_exp(), complex_array_sqrt();
  1178.   void complex_array_sin(), complex_array_cos();
  1179.   void complex_array_asin(), complex_array_acos(), complex_array_atan();
  1180.   PRIMITIVE_HEADER (3);
  1181.   CHECK_ARG (1, FIXNUM_P);    /* operation opcode */
  1182.   CHECK_ARG (2, ARRAY_P);    /* input array -- n      real part         */
  1183.   CHECK_ARG (3, ARRAY_P);    /* input array -- n      imag part         */
  1184.   n = ARRAY_LENGTH(ARG_REF(2));
  1185.   if (n != ARRAY_LENGTH(ARG_REF(3))) error_bad_range_arg(3);
  1186.   a  = ARRAY_CONTENTS(ARG_REF(2)); /*  real part */
  1187.   b  = ARRAY_CONTENTS(ARG_REF(3)); /*  imag part */
  1188.   opcode = arg_nonnegative_integer(1);
  1189.   if (opcode==1)
  1190.     complex_array_to_polar(a,b,n);
  1191.   else if (opcode==2)
  1192.     complex_array_exp(a,b,n);
  1193.   else if (opcode==3)
  1194.     complex_array_sqrt(a,b,n);
  1195.  
  1196.   else if (opcode==4)
  1197.     complex_array_sin(a,b,n);
  1198.   else if (opcode==5)
  1199.     complex_array_cos(a,b,n);
  1200.   /* for tan(z) use sin(z)/cos(z) */
  1201.   
  1202.   else if (opcode==6)
  1203.     complex_array_asin(a,b,n);
  1204.   else if (opcode==7)
  1205.     complex_array_acos(a,b,n);
  1206.   else if (opcode==8)
  1207.     complex_array_atan(a,b,n);
  1208.   
  1209.   else
  1210.     error_bad_range_arg(1);    /* illegal opcode */
  1211.   PRIMITIVE_RETURN (UNSPECIFIC);
  1212. }
  1213.  
  1214. void
  1215. complex_array_to_polar (a,b,n)
  1216.      REAL *a,*b;
  1217.      long n;
  1218. {
  1219.   long i;
  1220.   double x,y, temp;
  1221.   for (i=0; i<n; i++)
  1222.     {
  1223.       x = (double) a[i];
  1224.       y = (double) b[i];
  1225.       temp = sqrt(x*x + y*y);
  1226.       if (temp == 0.0)
  1227.     b[i] = 0.0;        /* choose angle = 0    for x,y=0,0 */
  1228.       else
  1229.     b[i] = (REAL) atan2(y,x);
  1230.       a[i]   = (REAL) temp;
  1231.     }
  1232. }
  1233.  
  1234. void
  1235. complex_array_exp (a,b,n)
  1236.      REAL *a,*b;
  1237.      long n;
  1238. {
  1239.   long i;
  1240.   double x,y, temp;
  1241.  
  1242.   for (i=0; i<n; i++)
  1243.     {
  1244.       x = (double) a[i];
  1245.       y = (double) b[i];
  1246.       if ((temp = exp(x)) == HUGE) error_bad_range_arg(2); /* overflow */
  1247.       a[i] = (REAL) (temp*cos(y));
  1248.       b[i] = (REAL) (temp*sin(y));
  1249.     }
  1250. }
  1251.  
  1252. void
  1253. complex_array_sqrt (a,b,n)
  1254.      REAL *a,*b;
  1255.      long n;
  1256. {
  1257.   long i;
  1258.   double x,y, r;
  1259.  
  1260.   for (i=0; i<n; i++)
  1261.     {
  1262.       x = (double) a[i];
  1263.       y = (double) b[i];
  1264.       r = sqrt( x*x + y*y);
  1265.       a[i] = sqrt((r+x)/2.0);
  1266.       if (y>=0.0)
  1267.     b[i] =  sqrt((r-x)/2.0); /* choose principal root */
  1268.       else            /* see Abramowitz (p.17 3.7.27) */
  1269.     b[i] = -sqrt((r-x)/2.0);
  1270.     }
  1271. }
  1272.  
  1273. void
  1274. complex_array_sin (a,b,n)
  1275.      REAL *a,*b;
  1276.      long n;
  1277. {
  1278.   long i;
  1279.   double x, ey,fy;
  1280.   REAL temp;
  1281.  
  1282.   for (i=0; i<n; i++)
  1283.     {
  1284.       x = (double) a[i];
  1285.       ey = exp((double) b[i]);    /* radius should be small to avoid overflow */
  1286.       fy = 1.0/ey;        /* exp(-y) */
  1287.       /* expanded (e(iz)-e(-iz))*(-.5i) formula */
  1288.       temp = (REAL) (sin(x) * (ey + fy) * 0.5);
  1289.       /* see my notes in Abram.p.71 */
  1290.       b[i] = (REAL) (cos(x) * (ey - fy) * 0.5);
  1291.       a[i] = temp;
  1292.     }
  1293. }
  1294.  
  1295. void
  1296. complex_array_cos (a,b,n)
  1297.      REAL *a,*b;
  1298.      long n;
  1299. {
  1300.   long i;
  1301.   double x, ey,fy;
  1302.   REAL temp;
  1303.  
  1304.   for (i=0; i<n; i++)
  1305.     {
  1306.       x = (double) a[i];
  1307.       ey = exp((double) b[i]);    /* radius should be small to avoid overflow */
  1308.       fy = 1.0/ey;        /* exp(-y) */
  1309.       /* expanded (e(iz)+e(-iz))*.5 formula */
  1310.       temp = (REAL) (cos(x) * (ey + fy) * 0.5);
  1311.       b[i] = (REAL) (sin(x) * (fy - ey) * 0.5); /* see my notes in Abram.p.71*/
  1312.       a[i] = temp;
  1313.     }
  1314. }
  1315.  
  1316.  
  1317. void
  1318. complex_array_asin (a,b,n)
  1319.      REAL *a,*b;
  1320.      long n;
  1321. { /* logarithmic formula as in R3.99, about 21ops plus log,atan - see my notes */
  1322.   long i; 
  1323.   double oldx,oldy, x,y, real,imag, r;
  1324.   
  1325.   for (i=0; i<n; i++)
  1326.   {
  1327.     oldx = (double) a[i];
  1328.     oldy = (double) b[i];
  1329.     
  1330.     x = 1.0 - oldx*oldx + oldy*oldy; /* 1 - z*z */
  1331.     y = -2.0 * oldx * oldy;
  1332.     
  1333.     r = sqrt(x*x + y*y);    /* sqrt(1-z*z)  */
  1334.     real = sqrt((r+x)/2.0);
  1335.     if (y>=0.0)
  1336.       imag =  sqrt((r-x)/2.0);    /* choose principal root */
  1337.     else            /* see Abramowitz (p.17 3.7.27) */
  1338.       imag = -sqrt((r-x)/2.0);
  1339.     
  1340.     real = real - oldy;        /* i*z + sqrt(...) */
  1341.     imag = imag + oldx;
  1342.     
  1343.     b[i] = (REAL) (- log (sqrt (real*real + imag*imag))); /* -i*log(...) */
  1344.     a[i] = (REAL) atan2( imag, real); /* chosen angle is okay 
  1345.                      Also 0/0 doesnot occur */
  1346.   }
  1347. }
  1348.  
  1349. void
  1350. complex_array_acos (a,b,n)
  1351.      REAL *a,*b;
  1352.      long n;
  1353. {
  1354.   long i;
  1355.  
  1356.   complex_array_asin (a,b,n);
  1357.   
  1358.   for (i=0; i<n; i++)
  1359.     {
  1360.       a[i] = PI_OVER_2 - a[i];
  1361.       b[i] =           - b[i];
  1362.     }
  1363. }
  1364.   
  1365.  
  1366. void
  1367. complex_array_atan (a,b,n)
  1368.      REAL *a,*b;
  1369.      long n;
  1370. { /* logarithmic formula, expanded, simplified - see my notes */
  1371.   long i; 
  1372.   double x,y, xx, real,imag, d;
  1373.   
  1374.   for (i=0; i<n; i++)
  1375.   {
  1376.     x = (double) a[i];
  1377.     y = (double) b[i];
  1378.     
  1379.     xx = x*x;
  1380.     imag = 1.0 + y;        /* temp var */
  1381.     d  = xx + imag*imag;
  1382.     
  1383.     real = (1 - y*y - xx) / d;
  1384.     imag = (2.0 * x)  / d;
  1385.     
  1386.     b[i] = (REAL) ((log (sqrt (real*real + imag*imag))) / -2.0);
  1387.     a[i] = (atan2 (imag,real)) / 2.0;
  1388.   }
  1389. }
  1390.  
  1391.  
  1392. /* complex-array-operation-1b!
  1393.    groups together procedures that use 1 complex-array & 1 number
  1394.    and store the result in place (e.g. invert 1/x) */
  1395.  
  1396. DEFINE_PRIMITIVE ("COMPLEX-ARRAY-OPERATION-1B!", Prim_complex_array_operation_1b, 4,4, 0)
  1397. {
  1398.   long n, i, opcode;
  1399.   REAL *a, *b, inf;
  1400.   void complex_array_invert ();
  1401.   PRIMITIVE_HEADER (4);
  1402.   CHECK_ARG (1, FIXNUM_P);    /* operation opcode */
  1403.   CHECK_ARG (2, ARRAY_P);    /* input array -- n      real part         */
  1404.   CHECK_ARG (3, ARRAY_P);    /* input array -- n      imag part         */
  1405.   inf = (arg_real (4));        /* User-Provided Infinity */
  1406.   n = ARRAY_LENGTH(ARG_REF(2));
  1407.   if (n != ARRAY_LENGTH(ARG_REF(3))) error_bad_range_arg(3);
  1408.   a  = ARRAY_CONTENTS(ARG_REF(2)); /*  real part */
  1409.   b  = ARRAY_CONTENTS(ARG_REF(3)); /*  imag part */
  1410.   opcode = arg_nonnegative_integer(1);
  1411.   if (opcode==1)
  1412.     complex_array_invert(a,b, n, inf);  /* performs 1/x */
  1413.   else if (opcode==2)
  1414.     error_bad_range_arg(1);    /* illegal opcode */
  1415.   else
  1416.     error_bad_range_arg(1);    /* illegal opcode */
  1417.   PRIMITIVE_RETURN (UNSPECIFIC);
  1418. }
  1419.  
  1420. void
  1421. complex_array_invert (a,b, n, inf)
  1422.      REAL *a,*b, inf;
  1423.      long n;
  1424. {
  1425.   long i;
  1426.   double x,y, r;
  1427.   for (i=0; i<n; i++)
  1428.     {
  1429.       x = (double) a[i];
  1430.       y = (double) b[i];
  1431.       r = (x*x + y*y);
  1432.       if (r==0.0)
  1433.     {
  1434.       a[i] = inf;
  1435.       b[i] = inf;
  1436.     }
  1437.       else
  1438.     {
  1439.       a[i] = (REAL)  x/r;
  1440.       b[i] = (REAL) -y/r;
  1441.     }
  1442.     }
  1443. }
  1444.  
  1445. /* complex-array-operation-1a
  1446.    groups together procedures that use 1 complex-array
  1447.    and store result in a 3rd real array. */
  1448.  
  1449. DEFINE_PRIMITIVE ("COMPLEX-ARRAY-OPERATION-1A", Prim_complex_array_operation_1a, 4,4, 0)
  1450. {
  1451.   long n, i, opcode;
  1452.   REAL *a, *b, *c;
  1453.   void complex_array_magnitude(), complex_array_angle();
  1454.   PRIMITIVE_HEADER (4);
  1455.   CHECK_ARG (1, FIXNUM_P);    /* operation opcode */
  1456.   CHECK_ARG (2, ARRAY_P);    /* input array -- n      real part         */
  1457.   CHECK_ARG (3, ARRAY_P);    /* input array -- n      imag part         */
  1458.   CHECK_ARG (4, ARRAY_P);    /* output array -- n                       */
  1459.   n = ARRAY_LENGTH(ARG_REF(2));
  1460.   if (n != ARRAY_LENGTH(ARG_REF(3))) error_bad_range_arg(3);
  1461.   if (n != ARRAY_LENGTH(ARG_REF(4))) error_bad_range_arg(4);
  1462.   a  = ARRAY_CONTENTS(ARG_REF(2)); /*  real part */
  1463.   b  = ARRAY_CONTENTS(ARG_REF(3)); /*  imag part */
  1464.   c  = ARRAY_CONTENTS(ARG_REF(4)); /*  output    */
  1465.   opcode = arg_nonnegative_integer(1);
  1466.   if (opcode==1)
  1467.     complex_array_magnitude(a,b,c,n);
  1468.   else if (opcode==2)
  1469.     complex_array_angle(a,b,c,n);
  1470.   else
  1471.     error_bad_range_arg(1);    /* illegal opcode */
  1472.   PRIMITIVE_RETURN (UNSPECIFIC);
  1473. }
  1474.  
  1475. void
  1476. complex_array_magnitude (a,b,c,n)
  1477.      REAL *a,*b,*c;
  1478.      long n;
  1479. {
  1480.   long i;
  1481.   for (i=0; i<n; i++)
  1482.     c[i] = (REAL) sqrt( (double) a[i]*a[i] + b[i]*b[i] );
  1483. }
  1484.  
  1485. void
  1486. complex_array_angle (a,b,c,n)
  1487.      REAL *a,*b,*c;
  1488.      long n;
  1489. {
  1490.   long i;
  1491.   for (i=0; i<n; i++)
  1492.   {
  1493.     if ((a[i] == 0.0) && (b[i]==0.0))
  1494.       c[i] = 0.0;        /* choose angle=0 for point (0,0) */
  1495.     else
  1496.       c[i] = (REAL) atan2( (double) b[i], (double) a[i]);
  1497.     /* angle ==   -pi (exclusive) to +pi (inclusive) */
  1498.   }
  1499. }
  1500.  
  1501. DEFINE_PRIMITIVE ("CS-ARRAY-MAGNITUDE!", Prim_cs_array_magnitude, 1, 1, 0)
  1502. {
  1503.   long n, i;
  1504.   REAL *a;
  1505.   void cs_array_magnitude();
  1506.   PRIMITIVE_HEADER (1);
  1507.   CHECK_ARG (1, ARRAY_P);
  1508.   a = ARRAY_CONTENTS(ARG_REF(1)); /* input cs-array */
  1509.   n = ARRAY_LENGTH(ARG_REF(1));    /* becomes a standard array on return  */
  1510.   cs_array_magnitude(a,n);
  1511.   PRIMITIVE_RETURN (UNSPECIFIC);
  1512. }
  1513.  
  1514. /* result is a standard array (even signal, real data) */
  1515. void
  1516. cs_array_magnitude (a,n)
  1517.      REAL *a;
  1518.      long n;
  1519. {
  1520.   long i, n2, ni;
  1521.   n2 = n/2;            /* integer division truncates down */
  1522.   a[0]  = (REAL) fabs((double) a[0]); /*   imag=0 */
  1523.   if (2*n2 == n)        /* even length, n2 is only real */
  1524.     a[n2] = (REAL) fabs((double) a[n2]); /*  imag=0 */
  1525.   else
  1526.     /* odd length, make the loop include the n2 index */
  1527.     n2 = n2+1;
  1528.   for (i=1; i<n2; i++)
  1529.     {
  1530.       ni = n-i;
  1531.       a[i]   = (REAL)  sqrt( (double) a[i]*a[i] + (double) a[ni]*a[ni] );
  1532.       a[ni]  = a[i];        /* even signal */
  1533.     }
  1534. }
  1535.  
  1536. /* Rectangular and Polar
  1537.    A cs-array has even magnitude and odd angle (almost)
  1538.    hence a polar cs-array stores magnitude in the first half (real part)
  1539.    and angle in the second half (imag part)
  1540.    except for a[0] real-only and a[n2] (n even)
  1541.    The angle of a[0] is either 0 (pos. sign) or pi (neg. sign),
  1542.    but there is no place in an n-point cs-array to store this, so
  1543.    a[0] and a[n2] when n even are left unchanged  when going polar.
  1544.    as opposed to taking their absolute values, for magnitude. */
  1545.  
  1546. DEFINE_PRIMITIVE ("CS-ARRAY-TO-POLAR!", Prim_cs_array_to_polar, 1,1, 0)
  1547. {
  1548.   long n, i;
  1549.   REAL *a;
  1550.   void cs_array_to_polar();
  1551.   PRIMITIVE_HEADER (1);
  1552.   CHECK_ARG (1, ARRAY_P);
  1553.   a  = ARRAY_CONTENTS(ARG_REF(1));
  1554.   n = ARRAY_LENGTH(ARG_REF(1));
  1555.   cs_array_to_polar(a,n);
  1556.   PRIMITIVE_RETURN (UNSPECIFIC);
  1557. }
  1558.  
  1559. void
  1560. cs_array_to_polar (a,n)
  1561.      REAL *a;
  1562.      long n;
  1563. {
  1564.   long i, n2;
  1565.   double real, imag;        /* temporary variables */
  1566.   n2 = n/2;            /* integer division truncates down */
  1567.   /* a[0] stores both magnitude and angle
  1568.      (pos. sign angle=0 , neg. sign angle=pi) */
  1569.   if (2*n2 == n)        /* even length, n2 is only real */
  1570.     ;                /* a[n2] stores sign information like a[0] */
  1571.   else
  1572.     /* odd length, make the loop include the n2 index */
  1573.     n2 = n2+1;
  1574.   for (i=1; i<n2; i++)
  1575.     {
  1576.       real = (double) a[i];
  1577.       imag = (double) a[n-i];
  1578.       a[i]   = (REAL)  sqrt( real*real + imag*imag );
  1579.       if (a[i] == 0.0)
  1580.     a[n-i] = 0.0;
  1581.       else
  1582.     a[n-i] = (REAL) atan2( imag, real );
  1583.     }
  1584. }
  1585.  
  1586. DEFINE_PRIMITIVE ("CS-ARRAY-TO-RECTANGULAR!", Prim_cs_array_to_rectangular, 1,1, 0)
  1587. {
  1588.   long n,n2, i;
  1589.   double magn,angl;        /* temporary variables */
  1590.   REAL *a;
  1591.   PRIMITIVE_HEADER (1);
  1592.   CHECK_ARG (1, ARRAY_P);
  1593.   a  = ARRAY_CONTENTS(ARG_REF(1));
  1594.   n = ARRAY_LENGTH(ARG_REF(1));
  1595.   n2 = n/2;            /* integer division truncates down */
  1596.   ;                /* a[0] is okay */
  1597.   if (2*n2 == n)        /* even length, n2 is real only */
  1598.     ;                /* a[n2] is okay */
  1599.   else
  1600.     n2 = n2+1; /* odd length, make the loop include the n2 index */
  1601.   for (i=1; i<n2; i++)
  1602.   { magn = (double) a[i];
  1603.     angl = (double) a[n-i];
  1604.     a[i]   = (REAL)  magn * cos(angl);
  1605.     a[n-i] = (REAL)  magn * sin(angl); }
  1606.   PRIMITIVE_RETURN (UNSPECIFIC);
  1607. }
  1608.  
  1609. /* Convolution in the Time-Domain */
  1610. /* In the following macro
  1611.    To1 and To2 should be (Length1-1) and (Length2-1) respectively. */
  1612.  
  1613. #define C_Convolution_Point_Macro(X, Y, To1, To2, N, Result)        \
  1614. {                                    \
  1615.   long Min_of_N_To1 = (min ((N), (To1)));                \
  1616.   long mi, N_minus_mi;                            \
  1617.   REAL Sum = 0.0;                            \
  1618.   for (mi = (max (0, ((N) - (To2)))), N_minus_mi = ((N) - mi);        \
  1619.        (mi <= Min_of_N_To1);                        \
  1620.        mi += 1, N_minus_mi -= 1)                    \
  1621.     Sum += ((X [mi]) * (Y [N_minus_mi]));                \
  1622.   (Result) = Sum;                            \
  1623. }
  1624.  
  1625. DEFINE_PRIMITIVE ("CONVOLUTION-POINT", Prim_convolution_point, 3, 3, 0)
  1626. {
  1627.   long Length1, Length2, N;
  1628.   REAL *Array1, *Array2;
  1629.   REAL C_Result;
  1630.   PRIMITIVE_HEADER (3);
  1631.   CHECK_ARG (1, ARRAY_P);
  1632.   CHECK_ARG (2, ARRAY_P);
  1633.   Length1 = ARRAY_LENGTH (ARG_REF (1));
  1634.   Length2 = ARRAY_LENGTH (ARG_REF (2));
  1635.   N = (arg_nonnegative_integer (3));
  1636.   Array1 = ARRAY_CONTENTS (ARG_REF (1));
  1637.   Array2 = ARRAY_CONTENTS (ARG_REF (2));
  1638.   C_Convolution_Point_Macro(Array1, Array2, Length1-1, Length2-1, N, C_Result);
  1639.   PRIMITIVE_RETURN (double_to_flonum ((double) C_Result));
  1640. }
  1641.  
  1642. DEFINE_PRIMITIVE ("ARRAY-CONVOLUTION-IN-TIME!", Prim_array_convolution_in_time, 3, 3, 0)
  1643. {
  1644.   long n,m,l, n_1,m_1, i;
  1645.   REAL *a,*b,*c;
  1646.   PRIMITIVE_HEADER (3);
  1647.   CHECK_ARG (1, ARRAY_P);
  1648.   CHECK_ARG (2, ARRAY_P);
  1649.   CHECK_ARG (3, ARRAY_P);
  1650.   a = ARRAY_CONTENTS(ARG_REF(1));
  1651.   b = ARRAY_CONTENTS(ARG_REF(2));
  1652.   c = ARRAY_CONTENTS(ARG_REF(3));
  1653.   n = ARRAY_LENGTH(ARG_REF(1));
  1654.   m = ARRAY_LENGTH(ARG_REF(2));
  1655.   l = n+m-1;            /* resulting length */
  1656.   if (l != ARRAY_LENGTH(ARG_REF(3))) error_bad_range_arg(3);
  1657.   n_1 = n-1; m_1 = m-1;
  1658.   for (i=0; i<l; i++)
  1659.   { C_Convolution_Point_Macro(a, b, n_1, m_1, i, c[i]); }
  1660.   PRIMITIVE_RETURN (UNSPECIFIC);
  1661. }
  1662.  
  1663. DEFINE_PRIMITIVE ("ARRAY-MULTIPLY-INTO-SECOND-ONE!", Prim_array_multiply_into_second_one, 2, 2, 0)
  1664. {
  1665.   long Length, i;
  1666.   REAL *To_Here;
  1667.   REAL *From_Here_1, *From_Here_2;
  1668.   PRIMITIVE_HEADER (2);
  1669.   CHECK_ARG (1, ARRAY_P);
  1670.   CHECK_ARG (2, ARRAY_P);
  1671.   Length = ARRAY_LENGTH (ARG_REF (1));
  1672.   if (Length != ARRAY_LENGTH (ARG_REF (2))) error_bad_range_arg (2);
  1673.   From_Here_1 = ARRAY_CONTENTS (ARG_REF (1));
  1674.   From_Here_2 = ARRAY_CONTENTS (ARG_REF (2));
  1675.   To_Here = From_Here_2;
  1676.   for (i=0; i < Length; i++)
  1677.     {
  1678.       *To_Here++ = (*From_Here_1) * (*From_Here_2);
  1679.       From_Here_1++ ;
  1680.       From_Here_2++ ;
  1681.     }
  1682.   PRIMITIVE_RETURN (UNSPECIFIC);
  1683. }
  1684.  
  1685. /* complex-array-operation-2!
  1686.    groups together procedures that use 2 complex-arrays
  1687.    and store result in either 1st or 2nd */
  1688.  
  1689. DEFINE_PRIMITIVE ("COMPLEX-ARRAY-OPERATION-2!", Prim_complex_array_operation_2, 5,5, 0)
  1690. {
  1691.   long n, opcode;
  1692.   REAL *ax,*ay, *bx,*by;
  1693.   void complex_array_multiply_into_second_one();
  1694.   PRIMITIVE_HEADER (5);
  1695.   CHECK_ARG (1, FIXNUM_P);    /* operation opcode */
  1696.   CHECK_ARG (2, ARRAY_P);    /* ax array -- n      real         */
  1697.   CHECK_ARG (3, ARRAY_P);    /* ay array -- n      imag         */
  1698.   CHECK_ARG (4, ARRAY_P);    /* bx array -- n      real         */
  1699.   CHECK_ARG (5, ARRAY_P);    /* by array -- n      imag         */
  1700.   n = ARRAY_LENGTH(ARG_REF(2));
  1701.   if (n != ARRAY_LENGTH(ARG_REF(3))) error_bad_range_arg(3);
  1702.   if (n != ARRAY_LENGTH(ARG_REF(4))) error_bad_range_arg(4);
  1703.   if (n != ARRAY_LENGTH(ARG_REF(4))) error_bad_range_arg(5);
  1704.   ax  = ARRAY_CONTENTS(ARG_REF(2)); /*  real */
  1705.   ay  = ARRAY_CONTENTS(ARG_REF(3)); /*  imag */
  1706.   bx  = ARRAY_CONTENTS(ARG_REF(4)); /*  real */
  1707.   by  = ARRAY_CONTENTS(ARG_REF(5)); /*  imag */
  1708.   opcode = arg_nonnegative_integer(1);
  1709.   if (opcode==1)
  1710.     complex_array_multiply_into_second_one(ax,ay,bx,by, n);
  1711.   else if (opcode==2)
  1712.     error_bad_range_arg(1);    /* illegal opcode */
  1713.   else
  1714.     error_bad_range_arg(1);    /* illegal opcode */
  1715.   PRIMITIVE_RETURN (UNSPECIFIC);
  1716. }
  1717.  
  1718. void
  1719. complex_array_multiply_into_second_one (ax,ay,bx,by, n)
  1720.      REAL *ax,*ay,*bx,*by;
  1721.      long n;
  1722. {
  1723.   long i;
  1724.   REAL temp;
  1725.   for (i=0;i<n;i++)
  1726.     {
  1727.       temp   = ax[i]*bx[i]  -  ay[i]*by[i]; /*  real part */
  1728.       by[i]  = ax[i]*by[i]  +  ay[i]*bx[i]; /*  imag part */
  1729.       bx[i]  = temp;
  1730.     }
  1731. }
  1732.  
  1733. void
  1734. C_Array_Complex_Multiply_Into_First_One (a,b,c,d, length) /* used in fft.c */
  1735.      REAL *a,*b,*c,*d;
  1736.      long length;
  1737. {
  1738.   long i;
  1739.   REAL temp;
  1740.   for (i=0;i<length;i++)
  1741.     {
  1742.       temp = a[i]*c[i] - b[i]*d[i];
  1743.       b[i] = a[i]*d[i] + b[i]*c[i];
  1744.       a[i] = temp;
  1745.     }
  1746. }
  1747.  
  1748. DEFINE_PRIMITIVE ("ARRAY-DIVIDE-INTO-XXX!", Prim_array_divide_into_xxx, 4,4, 0)
  1749. {
  1750.   long n, i, one_or_two;
  1751.   REAL *x,*y,*z, inf;
  1752.   void array_divide_into_z();
  1753.   PRIMITIVE_HEADER (4);
  1754.   CHECK_ARG (1, ARRAY_P);
  1755.   CHECK_ARG (2, ARRAY_P);
  1756.   inf = (arg_real (3));
  1757.   /* where to store result of division */
  1758.   one_or_two = (arg_nonnegative_integer (4));
  1759.   x = ARRAY_CONTENTS(ARG_REF(1));
  1760.   y = ARRAY_CONTENTS(ARG_REF(2));
  1761.   n = ARRAY_LENGTH(ARG_REF(1));
  1762.   if (n!=(ARRAY_LENGTH(ARG_REF(2)))) error_bad_range_arg(2);
  1763.   if (one_or_two == 1)
  1764.     array_divide_into_z( x,y, x,  n, inf);
  1765.   else if (one_or_two == 2)
  1766.     array_divide_into_z( x,y, y,  n, inf);
  1767.   else
  1768.     error_bad_range_arg(4);
  1769.   PRIMITIVE_RETURN (UNSPECIFIC);
  1770. }
  1771.  
  1772. void
  1773. array_divide_into_z (x,y, z, n, inf) /* z can either x or y */
  1774.      REAL *x,*y,*z, inf;
  1775.      long n;
  1776. {
  1777.   long i;
  1778.   for (i=0; i<n; i++)
  1779.     {
  1780.       if (y[i] == 0.0)
  1781.     {
  1782.       if (x[i] == 0.0)
  1783.         z[i] = 1.0;
  1784.       else
  1785.         z[i] = inf * x[i];
  1786.     }
  1787.       else
  1788.     z[i] = x[i] / y[i];
  1789.     }
  1790. }
  1791.  
  1792. /* complex-array-operation-2b!
  1793.    groups together procedures that use 2 complex-arrays
  1794.    & 1 additional real number
  1795.    and store result in either 1st or 2nd (e.g. division) */
  1796.  
  1797. DEFINE_PRIMITIVE ("COMPLEX-ARRAY-OPERATION-2B!", Prim_complex_array_operation_2b, 6,6, 0)
  1798. { long n, opcode;
  1799.   REAL *ax,*ay, *bx,*by,  inf;
  1800.   void complex_array_divide_into_z();
  1801.   PRIMITIVE_HEADER (6);
  1802.   CHECK_ARG (1, FIXNUM_P);    /* operation opcode */
  1803.   CHECK_ARG (2, ARRAY_P);    /* ax array -- n      real         */
  1804.   CHECK_ARG (3, ARRAY_P);    /* ay array -- n      imag         */
  1805.   CHECK_ARG (4, ARRAY_P);    /* bx array -- n      real         */
  1806.   CHECK_ARG (5, ARRAY_P);    /* by array -- n      imag         */
  1807.   inf = (arg_real (6));        /* User-Provided Infinity */
  1808.   n = ARRAY_LENGTH(ARG_REF(2));
  1809.   if (n != ARRAY_LENGTH(ARG_REF(3))) error_bad_range_arg(3);
  1810.   if (n != ARRAY_LENGTH(ARG_REF(4))) error_bad_range_arg(4);
  1811.   if (n != ARRAY_LENGTH(ARG_REF(4))) error_bad_range_arg(5);
  1812.   ax  = ARRAY_CONTENTS(ARG_REF(2)); /*  real */
  1813.   ay  = ARRAY_CONTENTS(ARG_REF(3)); /*  imag */
  1814.   bx  = ARRAY_CONTENTS(ARG_REF(4)); /*  real */
  1815.   by  = ARRAY_CONTENTS(ARG_REF(5)); /*  imag */
  1816.   opcode = arg_nonnegative_integer(1);
  1817.   if (opcode==1)
  1818.     complex_array_divide_into_z (ax,ay,bx,by, ax,ay,  n, inf);
  1819.   else if (opcode==2)
  1820.     complex_array_divide_into_z (ax,ay,bx,by, bx,by,  n, inf);
  1821.   else
  1822.     error_bad_range_arg(1);    /* illegal opcode */
  1823.   PRIMITIVE_RETURN (UNSPECIFIC);
  1824. }
  1825.  
  1826. void
  1827. complex_array_divide_into_z (xr,xi, yr,yi, zr,zi, n, inf)
  1828.      REAL *xr,*xi, *yr,*yi, *zr,*zi, inf;
  1829.      long n;
  1830. {
  1831.   long i;
  1832.   fast double temp, radius;
  1833.   for (i=0; i<n; i++)
  1834.     {
  1835.       radius = (double) (yr[i] * yr[i]) + (yi[i] * yi[i]); /* denominator */
  1836.       if (radius == 0.0)
  1837.     {
  1838.       if (xr[i] == 0.0)
  1839.         zr[i] = 1.0;
  1840.       else
  1841.         zr[i] = inf * xr[i];
  1842.       if (xi[i] == 0.0)
  1843.         zi[i] = 1.0;
  1844.       else
  1845.         zi[i] = inf * xi[i];
  1846.     }
  1847.       else
  1848.     {
  1849.       temp =  (double) (xr[i] * yr[i]  +  xi[i] * yi[i]);
  1850.       zi[i] = (REAL) (xi[i] * yr[i]  -  xr[i] * yi[i]) / radius;
  1851.       zr[i] = (REAL) temp                              / radius;
  1852.     }
  1853.     }
  1854. }
  1855.  
  1856. DEFINE_PRIMITIVE ("ARRAY-LINEAR-SUPERPOSITION-INTO-SECOND-ONE!", Prim_array_linear_superposition_into_second_one, 4, 4, 0)
  1857. {
  1858.   long n, i;
  1859.   REAL *To_Here, Coeff1, Coeff2;
  1860.   REAL *From_Here_1, *From_Here_2;
  1861.   PRIMITIVE_HEADER (4);
  1862.   Coeff1 = (arg_real (1));
  1863.   CHECK_ARG (2, ARRAY_P);
  1864.   Coeff2 = (arg_real (3));
  1865.   CHECK_ARG (4, ARRAY_P);
  1866.   n = (ARRAY_LENGTH (ARG_REF (2)));
  1867.   if (n != (ARRAY_LENGTH (ARG_REF (4))))
  1868.     error_bad_range_arg (4);
  1869.   From_Here_1 = (ARRAY_CONTENTS (ARG_REF (2)));
  1870.   From_Here_2 = (ARRAY_CONTENTS (ARG_REF (4)));
  1871.   To_Here = From_Here_2;
  1872.   for (i=0; i < n; i++)
  1873.     {
  1874.       *To_Here++ = (Coeff1 * (*From_Here_1)) + (Coeff2 * (*From_Here_2));
  1875.       From_Here_1++ ;
  1876.       From_Here_2++ ;
  1877.     }
  1878.   PRIMITIVE_RETURN (UNSPECIFIC);
  1879. }
  1880.  
  1881. /*  m_pi = 3.14159265358979323846264338327950288419716939937510; */
  1882.  
  1883. DEFINE_PRIMITIVE ("SAMPLE-PERIODIC-FUNCTION", Prim_sample_periodic_function, 4, 4, 0)
  1884. {
  1885.   long N, i, Function_Number;
  1886.   double Signal_Frequency, Sampling_Frequency, DT, DTi;
  1887.   double twopi = 6.28318530717958;
  1888.   SCHEME_OBJECT Result, Pfunction_number, Psignal_frequency;
  1889.   SCHEME_OBJECT Pfunction_Number;
  1890.   REAL *To_Here;
  1891.   double unit_square_wave(), unit_triangle_wave();
  1892.   PRIMITIVE_HEADER (4);
  1893.   Function_Number = (arg_index_integer (1, 11));
  1894.   Signal_Frequency = (arg_real_number (2));
  1895.   if (Signal_Frequency == 0)
  1896.     error_bad_range_arg (2);
  1897.   Sampling_Frequency = (arg_real_number (3));
  1898.   if (Sampling_Frequency == 0)
  1899.     error_bad_range_arg (3);
  1900.   N = (arg_nonnegative_integer (4));
  1901.   Result = (allocate_array (N));
  1902.   To_Here = ARRAY_CONTENTS(Result);
  1903.   DT = (double) (twopi * Signal_Frequency * (1 / Sampling_Frequency));
  1904.   if (Function_Number == 0)
  1905.     for (i=0, DTi=0.0; i < N; i++, DTi += DT)
  1906.       *To_Here++ = (REAL) cos(DTi);
  1907.   else if (Function_Number == 1)
  1908.     for (i=0, DTi=0.0; i < N; i++, DTi += DT)
  1909.       *To_Here++ = (REAL) sin(DTi);
  1910.   else if (Function_Number == 2)
  1911.     for (i=0, DTi=0.0; i < N; i++, DTi += DT)
  1912.       *To_Here++ = (REAL) unit_square_wave(DTi);
  1913.   else if (Function_Number == 3)
  1914.     for (i=0, DTi=0.0; i < N; i++, DTi += DT)
  1915.       *To_Here++ = (REAL) unit_triangle_wave(DTi);
  1916.   else
  1917.     error_bad_range_arg (1);
  1918.   PRIMITIVE_RETURN (Result);
  1919. }
  1920. ;
  1921. double
  1922. hamming (t, length)
  1923.      double t, length;
  1924. {
  1925.   double twopi = 6.28318530717958;
  1926.   double pi = twopi/2.;
  1927.   double t_bar = cos(twopi * (t / length));
  1928.   if ((t<length) && (t>0.0)) return(.08 + .46 * (1 - t_bar));
  1929.   else return (0);
  1930. }
  1931.  
  1932. double
  1933. unit_square_wave (t)
  1934.      double t;
  1935. {
  1936.   double twopi = 6.28318530717958;
  1937.   double fmod(), fabs();
  1938.   double pi = twopi/2.;
  1939.   double t_bar = ((REAL) fabs(fmod( ((double) t), twopi)));
  1940.   if (t_bar < pi) return(1);
  1941.   else return(-1);
  1942. }
  1943.  
  1944. double
  1945. unit_triangle_wave (t)
  1946.      double t;
  1947. {
  1948.   double twopi = 6.28318530717958;
  1949.   double pi = twopi/2.;
  1950.   double pi_half = pi/2.;
  1951.   double three_pi_half = pi+pi_half;
  1952.   double t_bar = ((double) fabs(fmod( ((double) t), twopi)));
  1953.   if (t_bar<pi_half)
  1954.     return (-(t_bar/pi));
  1955.   else if (t_bar<pi)
  1956.     return (t_bar/pi);
  1957.   else if (t_bar<three_pi_half)
  1958.     return ((twopi-t_bar)/pi);
  1959.   else
  1960.     return (-((twopi-t_bar)/pi));
  1961. }
  1962.  
  1963. DEFINE_PRIMITIVE ("ARRAY-HANNING!", Prim_array_hanning, 2,2, 0)
  1964. {
  1965.   long n, hanning_power;
  1966.   REAL *a;
  1967.   void C_Array_Make_Hanning();
  1968.   PRIMITIVE_HEADER (2);
  1969.   CHECK_ARG (1, ARRAY_P);    /* input array -- n */
  1970.   CHECK_ARG (2, FIXNUM_P);    /* hanning power */
  1971.   a  = ARRAY_CONTENTS(ARG_REF(1));
  1972.   n = ARRAY_LENGTH(ARG_REF(1));
  1973.   hanning_power = arg_nonnegative_integer(2);
  1974.   C_Array_Make_Hanning (a, n, hanning_power);
  1975.   PRIMITIVE_RETURN (UNSPECIFIC);
  1976. }
  1977.  
  1978. void
  1979. C_Array_Make_Hanning (f1, length, power)
  1980.      REAL f1[];
  1981.      long length, power;
  1982. {
  1983.   double window_length;
  1984.   long i;
  1985.   double integer_power(), hanning();
  1986.   window_length = ((double) length);
  1987.   for (i=0;i<length;i++)
  1988.     {
  1989.       f1[i] = ((REAL) hanning(((double) i), window_length));
  1990.       f1[i] = (REAL) integer_power(((double) f1[i]), power);
  1991.     }
  1992. }
  1993.  
  1994. double
  1995. hanning (t, length)
  1996.      double t, length;
  1997. {
  1998.   double twopi = 6.283185307179586476925287;
  1999.   double t_bar;
  2000.   t_bar = cos(twopi * (t / length));
  2001.   if ((t<length) && (t>0.0))
  2002.     return(.5 * (1 - t_bar));
  2003.   else
  2004.     return (0.0);
  2005. }
  2006.  
  2007. double
  2008. integer_power (a, n)
  2009.      double a;
  2010.      long n;
  2011. {
  2012.   double b;
  2013.   double integer_power();
  2014.  
  2015.   if (n<0) exit(-1);
  2016.   else if (n==0) return(1.0);
  2017.   else if (n==1) return(a);
  2018.   else if ((n%2) == 0)
  2019.     {
  2020.       b = integer_power(a, n/2);
  2021.       return(b*b);
  2022.     }
  2023.   else
  2024.     return(a * integer_power(a, (n-1)));
  2025. }
  2026.  
  2027. /* array-operation-1!
  2028.    groups together procedures that use 1 array
  2029.    and store the result in place (e.g. random) */
  2030.  
  2031. DEFINE_PRIMITIVE ("ARRAY-OPERATION-1!", Prim_array_operation_1, 2,2, 0)
  2032. {
  2033.   long n, opcode;
  2034.   REAL *a;
  2035.   void array_random();
  2036.   PRIMITIVE_HEADER (2);
  2037.   CHECK_ARG (1, FIXNUM_P);    /* operation opcode */
  2038.   CHECK_ARG (2, ARRAY_P);    /* input array -- n */
  2039.   n = ARRAY_LENGTH(ARG_REF(2));
  2040.   a  = ARRAY_CONTENTS(ARG_REF(2));
  2041.   opcode = arg_nonnegative_integer(1);
  2042.   if (opcode==1)
  2043.     array_random(a,n);
  2044.   else if (opcode==2)
  2045.     error_bad_range_arg(1);    /* illegal opcode */
  2046.   else
  2047.     error_bad_range_arg(1);    /* illegal opcode */
  2048.   PRIMITIVE_RETURN (UNSPECIFIC);
  2049. }
  2050.  
  2051. void
  2052. array_random (a,n)
  2053.      REAL *a;
  2054.      long n;
  2055. {
  2056.   long i;
  2057.   /* HPUX 3: Rand uses a multiplicative congruential random-number generator
  2058.      with period 2^32 that returns successive pseudo-random numbers in the
  2059.      range from 0 to 2^15-1 */
  2060.   for (i=0;i<n;i++)
  2061.     a[i] = ((REAL) rand()) * (3.0517578125e-5);
  2062.   /* 3.051xxx = 2^(-15) makes the range from 0 to 1 */
  2063. }
  2064.  
  2065. /* The following should go away.
  2066.    superceded by ARRAY-CONS-INTEGERS, ARRAY-UNARY-FUNCTION and array-random */
  2067. DEFINE_PRIMITIVE ("SAMPLE-APERIODIC-FUNCTION", Prim_sample_aperiodic_function, 3, 3, 0)
  2068. {
  2069.   long N, i, Function_Number;
  2070.   double Sampling_Frequency, DT, DTi;
  2071.   double twopi = 6.28318530717958;
  2072.   SCHEME_OBJECT Result;
  2073.   REAL *To_Here, twopi_dt;
  2074.   PRIMITIVE_HEADER (3);
  2075.   Function_Number = (arg_index_integer (1, 7));
  2076.   Sampling_Frequency = (arg_real_number (2));
  2077.   if (Sampling_Frequency == 0)
  2078.     error_bad_range_arg (2);
  2079.   N = (arg_nonnegative_integer (3));
  2080.   Result = (allocate_array (N));
  2081.   To_Here = ARRAY_CONTENTS(Result);
  2082.   DT = (twopi * (1 / Sampling_Frequency));
  2083.   if      (Function_Number == 0)
  2084.     /* HPUX 3: Rand uses a multiplicative congruential random-number generator
  2085.        with period 2^32 that returns successive pseudo-random numbers in the
  2086.        range from 0 to 2^15-1 */
  2087.     for (i=0; i<N; i++)
  2088.       /* 2^(-15) makes range from 0 to 1 */
  2089.       *To_Here++ = 3.0517578125e-5 * ((REAL) rand());
  2090.   else if (Function_Number == 1)
  2091.   { double length=DT*N;
  2092.     for (i=0, DTi=0.0; i < N; i++, DTi += DT)
  2093.       *To_Here++ = (REAL) hanning(DTi, length);
  2094.   }
  2095.   else if (Function_Number == 2)
  2096.   { double length=DT*N;
  2097.     for (i=0, DTi=0.0; i < N; i++, DTi += DT)
  2098.       *To_Here++ = (REAL) hamming(DTi, length);
  2099.   }
  2100.   else if (Function_Number == 3)
  2101.     for (i=0, DTi=0.0; i < N; i++, DTi += DT)
  2102.       *To_Here++ = (REAL) sqrt(DTi);
  2103.   else if (Function_Number == 4)
  2104.     for (i=0, DTi=0.0; i < N; i++, DTi += DT)
  2105.       *To_Here++ = (REAL) log(DTi);
  2106.   else if (Function_Number == 5)
  2107.     for (i=0, DTi=0.0; i < N; i++, DTi += DT)
  2108.       *To_Here++ = (REAL) exp(DTi);
  2109.   else
  2110.     error_bad_range_arg (1);
  2111.   PRIMITIVE_RETURN (Result);
  2112. }
  2113.  
  2114. DEFINE_PRIMITIVE ("ARRAY-PERIODIC-DOWNSAMPLE", Prim_array_periodic_downsample, 2, 2, 0)
  2115. {
  2116.   long Length, Pseudo_Length, Sampling_Ratio;
  2117.   REAL *Array, *To_Here;
  2118.   SCHEME_OBJECT Result;
  2119.   long i, array_index;
  2120.   PRIMITIVE_HEADER (2);
  2121.   CHECK_ARG (1, ARRAY_P);
  2122.   Length = ARRAY_LENGTH(ARG_REF (1));
  2123.   Sampling_Ratio = ((arg_integer (2)) % Length);
  2124.   if (Sampling_Ratio < 1)
  2125.     error_bad_range_arg (2);
  2126.   Array = ARRAY_CONTENTS(ARG_REF (1));
  2127.   Result = (allocate_array (Length));
  2128.   To_Here = ARRAY_CONTENTS(Result);
  2129.   Pseudo_Length = Length * Sampling_Ratio;
  2130.   /* new Array has the same Length by assuming periodicity */
  2131.   for (i=0; i<Pseudo_Length; i += Sampling_Ratio)
  2132.   { array_index = i % Length;
  2133.     *To_Here++ = Array[array_index]; }
  2134.   PRIMITIVE_RETURN (Result);
  2135. }
  2136.  
  2137. /* Shift is not done in place (no side-effects). */
  2138. DEFINE_PRIMITIVE ("ARRAY-PERIODIC-SHIFT", Prim_array_periodic_shift, 2, 2, 0)
  2139. {
  2140.   long Length, Shift;
  2141.   REAL *Array, *To_Here;
  2142.   SCHEME_OBJECT Result;
  2143.   long i, array_index;
  2144.   PRIMITIVE_HEADER (2);
  2145.   CHECK_ARG (1, ARRAY_P);
  2146.   Length = ARRAY_LENGTH(ARG_REF (1));
  2147.   /* periodic waveform, same sign as dividend */
  2148.   Shift = ((arg_integer (2)) % Length);
  2149.   Array = ARRAY_CONTENTS(ARG_REF (1));
  2150.   Result = (allocate_array (Length));
  2151.   To_Here = ARRAY_CONTENTS(Result);
  2152.   /* new Array has the same Length by assuming periodicity */
  2153.   for (i=0; i<Length; i++) {
  2154.     array_index = (i+Shift) % Length;
  2155.     if (array_index<0) array_index = Length + array_index; /* wrap around */
  2156.     *To_Here++ = Array[array_index]; }
  2157.   PRIMITIVE_RETURN (Result);
  2158. }
  2159.  
  2160. /* This is done here because array-map is very slow */
  2161. DEFINE_PRIMITIVE ("ARRAY-APERIODIC-DOWNSAMPLE", Prim_array_aperiodic_downsample, 2, 2, 0)
  2162. {
  2163.   long Length, New_Length, Sampling_Ratio;
  2164.   REAL *Array, *To_Here;
  2165.   SCHEME_OBJECT Result;
  2166.   long i, array_index;
  2167.   PRIMITIVE_HEADER (2);
  2168.   CHECK_ARG (1, ARRAY_P);
  2169.   Array = ARRAY_CONTENTS(ARG_REF (1));
  2170.   Length = ARRAY_LENGTH(ARG_REF (1));
  2171.   Sampling_Ratio = (arg_integer_in_range (2, 1, (Length + 1)));
  2172.   if (Length < 1) error_bad_range_arg (1);
  2173.   /* 1 for first one and then the rest --- integer division chops */
  2174.   New_Length = 1 + ((Length-1)/Sampling_Ratio);
  2175.   Result = (allocate_array (New_Length));
  2176.   To_Here = ARRAY_CONTENTS(Result);
  2177.   for (i=0; i<Length; i += Sampling_Ratio)
  2178.     *To_Here++ = Array[i];
  2179.   PRIMITIVE_RETURN (Result);
  2180. }
  2181.  
  2182. /* one more hack for speed */
  2183.  
  2184. /* (SOLVE-SYSTEM A B N)
  2185.     Solves the system of equations Ax = b.  A and B are
  2186.     arrays and b is the order of the system.  Returns x.
  2187.     From the Fortran procedure in Strang. */
  2188.  
  2189. DEFINE_PRIMITIVE ("SOLVE-SYSTEM", Prim_gaussian_elimination, 2, 2, 0)
  2190. {
  2191.   REAL *A, *B, *X;
  2192.   long Length;
  2193.   SCHEME_OBJECT Result;
  2194.   PRIMITIVE_HEADER (2);
  2195.   CHECK_ARG (1, ARRAY_P);
  2196.   CHECK_ARG (2, ARRAY_P);
  2197.   Length  = ARRAY_LENGTH(ARG_REF (2));
  2198.   if ((Length*Length) != ARRAY_LENGTH(ARG_REF (1))) error_bad_range_arg (2);
  2199.   A = ARRAY_CONTENTS(ARG_REF (1));
  2200.   B = ARRAY_CONTENTS(ARG_REF (2));
  2201.   Result = (allocate_array (Length));
  2202.   X = ARRAY_CONTENTS(Result);
  2203.   C_Array_Copy(B, X, Length);
  2204.   C_Gaussian_Elimination(A, X, Length);
  2205.   PRIMITIVE_RETURN (Result);
  2206. }
  2207.  
  2208. /* C routine side-effects b. */
  2209. C_Gaussian_Elimination (a, b, n)
  2210.      REAL *a, *b;
  2211.      long n;
  2212. {
  2213.   long *pvt;
  2214.   REAL p, t;
  2215.   long i, j, k, m;
  2216.   Primitive_GC_If_Needed (n);
  2217.   pvt = ((long *) Free);
  2218.   *(pvt+n-1) = 1;
  2219.   if (n != 1) {
  2220.     for (k=1; k<n; k++) {
  2221.       m = k;
  2222.       for (i=k+1; i<=n; i++)
  2223.     if (fabs(*(a+i+(k-1)*n-1)) > fabs(*(a+m+(k-1)*n-1)))
  2224.       m = i;
  2225.       *(pvt+k-1) = m;
  2226.       if (m != k)
  2227.     *(pvt+n-1) = - *(pvt+n-1);
  2228.       p = *(a+m+(k-1)*n-1);
  2229.       *(a+m+(k-1)*n-1) = *(a+k+(k-1)*n-1);
  2230.       *(a+k+(k-1)*n-1) = p;
  2231.       if (p != 0.0) {
  2232.     for (i=k+1; i<=n; i++)
  2233.       *(a+i+(k-1)*n-1) = - *(a+i+(k-1)*n-1) / p;
  2234.     for (j=k+1; j<=n; j++) {
  2235.       t = *(a+m+(j-1)*n-1);
  2236.       *(a+m+(j-1)*n-1) = *(a+k+(j-1)*n-1);
  2237.       *(a+k+(j-1)*n-1) = t;
  2238.       if (t != 0.0)
  2239.         for (i=k+1; i<=n; i++)
  2240.           *(a+i+(j-1)*n-1) = *(a+i+(j-1)*n-1) + *(a+i+(k-1)*n-1) * t;
  2241.     }
  2242.       }
  2243.     }
  2244.     for (k=1; k<n; k++) {
  2245.       m = *(pvt+k-1);
  2246.       t = *(b+m-1);
  2247.       *(b+m-1) = *(b+k-1);
  2248.       *(b+k-1) = t;
  2249.       for (i=k+1; i<=n; i++)
  2250.     *(b+i-1) = *(b+i-1) + *(a+i+(k-1)*n-1) * t;
  2251.     }
  2252.     for (j=1; j<n; j++) {
  2253.       k = n - j + 1;
  2254.       *(b+k-1) = *(b+k-1) / *(a+k+(k-1)*n-1);
  2255.       t = - *(b+k-1);
  2256.       for (i=1; i <= n-j; i++)
  2257.     *(b+i-1) = *(b+i-1) + *(a+i+(k-1)*n-1) * t;
  2258.     }
  2259.   }
  2260.   *b = *b / *a;
  2261.   return;
  2262. }
  2263.